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Four  topics  concerning  statistical  lediniques  for  the  analysis  of  multlresponse  ex- 
periments arc  presented.  These  techniques  can  be  used  to  (1)  set  up  a series  of 
multireaponse  experiments  (design)  that  will  yield  adequate  and  reliable  predictions 
of  the  responses  of  interest  on  the  basis  of  a certain  hypothesized  model;  (2)  fit  the 
hypothesized  multlresponse  model  using  data  collected  according  to  the  design  chosen 
in  (1);  and  (3]  determine  optimum  operating  conditions  on  the  input  variables  that 
lead  to  the  simultaneous  optima  of  the  predicted  values  of  the  responses,  within  a 
certain  region  of  interest  The  first  topic  is  an  extension  of  the  optimization  procedure 
discussed  by  Khuri  aud  Conlou  (1981).  The  simultaneous  optimization  of  multire- 
spouse  functions  represented  by  polynomial  regression  models  that  are  oot  neceasariiy 
of  the  same  form  and  degree  are  considered.  The  optimization  procedure  involves  the 
minimization  of  a distance  function  that  measures  the  distance  bctvreeu  the  vector  of 
estimated  responses  and  the  assumed  'ideal’  optimum  over  the  experimental  region. 
'The  second  topic  deals  with  multiresponse  models  with  block  effects  that  are  either 
fixed  or  random.  The  effect  of  blocking  oo  the  estimation  of  the  mean  response,  the 


of  the 


prediction  varinoce,  and  on  the  Bimultaneous  optima 
function  is  discussed.  The  final  topic  is  the  coiuilruction  of  designs  to  minimiae  a 
muitiresponse  version  of  the  integrated  mean-squared  error.  This  is  an  extension  of 
the  work  done  by  Kim  and  Draper  (1994)  to  cases  where  more  than  two  responses  and 

variance-covariance  matrix  associated  with  the  responses  is  unknown.  A sequetitiai 
procedure  wherein  the  design  pdms  arc  chosen  iterativeiy  is  presented. 


CHAPTER  1 


INTRODUCTION 

l.l  An  Overview  of  the  Preseni  Wnck 

Much  has  beec  writteD  concemiai  response  surface  methodology  (RSM)  involv- 
ing a single  response.  The  review  article  by  Myers  et  al.  (1989)  shows  how  much 
work  has  been  done  in  this  particular  area.  However,  the  same  cannot  be  said  about 
RSM  involving  several  responses.  The  design  and  analysis  of  experiments  with  mul- 

many  experiments  it  is  not  uncommon  to  measure  more  than  one  response  from  the 
same  experimental  unit.  These  experiments  are  called  mullirespoaie  espenments. 
For  example,  in  a chemical  process,  then  are  almost  always  several  properties  of  the 
product  that  are  of  iotccest,  such  as  the  yield  and  the  cost  of  making  the  product.  In 
the  development  of  food  products,  the  food  technologist  might  consider  in  addition 
to  the  Ravor  of  the  product,  the  shelf  life  and  subsequent  undesirable  byproducts.  In 
the  pharmaceutical  area,  a scientist  is  not  only  concerned  about  the  efficacy  of  the 
drug  but  also  of  the  possibility  of  serious  side  effects  ou  the  patient.  These  are  but 

advances  in  technology,  obtaining  data  that  describe  the  different  facets  of  a system 
or  product  is  becoming  mon  commonplace.  Thus  it  is  very  important  that  multire- 


It  is  therefore  the  purpose  of  this  dissertation  to  introduce  new  multirespoose  tech- 
niques which  pr^otly  are  missing  in  the  literature  on  the  design  and  analysis  of 
multirespoose  experiments. 

When  several  responses  are  considered  simultaneously,  they  should  not  be  treated 
as  Independent,  particularly  when  it  issuapecled  they  are  correlated.  In  multirespoose 
optimization,  for  example,  In  order  to  determine  optimum  operating  conditions,  it  is 
not  sufficient  to  only  determine  individual  optima  of  the  responses.  This  is  because 
optimal  conditions  for  one  response  may  be  far  from  optimal  or  even  physically  im- 
practical for  the  other  responses.  Also,  In  the  design  and  analysis  of  multirasponse 
experiments,  a design  that  is  satiafactory  for  one  response  may  produce  uosatiafactory 
results  for  the  rest  of  the  responses.  Thus  any  design  criterion  used  should  be  based 
on  the  consideration  of  all  Che  responses  together  rather  than  individually. 

One  important  objective  of  multiresponse  experimentation  is  tbe  determination 
of  optimum  operating  conditions  on  the  input  (control)  variables  that  will  lead  to  the 
simultaneous  optimization  of  tlie  predicted  values  of  the  responses  within  a certain 
region  of  interest.  For  example,  in  a tool  life  problem,  the  experimenter  may  wish  to 
determine  which  cutting  speed  and  depth  of  cut  will  maximize  the  lifetime  of  the  tool 
and  also  at  the  same  time  maximize  the  rate  of  material  removed.  Another  example 
is  the  chemical  process  situation  described  above.  The  scientist  may  want  to  find  the 
levels  of  temperature  and  pressure  which  butli  maxiuiize  the  yield  and  minimize  tbe 
cost  of  the  process.  There  are  many  examples  of  practical  situations  where  multiple 
responses  are  observed  and  where  there  is  a need  to  reach  some  type  of  compromise 
as  far  as  optimum  conditions  are  concerned.  Several  approaches  arc  available  in  the 
staCiatica  literature  which  can  be  used  to  determine  these  oplimaj  conditions.  Thsss 
approaches  include  the  ones  given  by  Myers  and  Carter  (1973),  Derrbiger  and  Suich 
(1980),  Khuri  and  Cunlon  (1981)  and  more  recently,  the  one  by  Del  Castillo  and 
Montgomery  (1993).  A review  of  these  procedures  is  given  in  Chapter  2.  In  Chapter 


is  introduced.  This  procedure  is 


3,  another  miilliresponse  optimisation  procedure 
an  exunslon  of  the  vwtrk  done  by  Khuri  and  Conlon  (1981)  to  the  case  where  the 
response  models  proposed  are  not  of  the  same  form  and  degree. 

There  are  many  experimental  situations  witere  it  is  impossible  to  perform  all  the 
runs  under  homogeneous  conditions.  Examples  of  homogeneous  conditions  might  be 
a single  period  of  time,  say  one  day,  a single  batch  of  homogeneous  taw  material,  or 
using  a single  machine  or  operator  In  a particular  process.  Multiple  days  or  batches  of 
raw  material  represent  possible  extraneous  sources  of  variation  (day,  batcli,  etc.)  and 
are  called  block  effects  Such  blocking  conditions  need  to  be  considered  in  the  analysis 
of  the  multireeponse  experiment.  Blocking  may  also  be  used  in  instances  where  it 
may  be  better  to  deliberately  vary  the  experimental  conditions  to  ensure  Chat  the 

practice.  For  example,  a chemical  engineer  may  run  a pilot  plant  experiment  with 

material  will  be  used  in  the  full-scale  process.  However,  very  little  infonnacion  is 
available  on  Che  analysis  of  multiresponse  experiments  with  block  effects.  Results 

I9€ta,b)  and  Khuri  (1993,1994,1996b)  but  nothing  much  has  been  done  about  the 
mullircsponse  case.  The  importance  of  considering  block  effects  in  single-response 
experimenta  ha.s  long  been  established  and  the  same  should  be  done  for  muitlresponse 
experiments.  A review  of  blocking  in  RSM  Is  given  in  Chapter  2.  Chapters  4 end  5 
deal  with  the  analysis  of  multiresponse  models  with  block  effects.  The  results  given  in 
these  two  chapters  are  the  multireeponse  generalizations  of  the  works  done  by  Khuri 
(1992,  1994,  1996b). 

Prior  to  running  a multiresponse  experiment,  the  researcher  must  hist  decide 
which  experimental  deeign  to  use.  The  task  of  selecting  a suitable  response  surface 


i easy  one.  There 


Dtal  designs  available  ic 


the  literature,  .^s  a result,  Myers  and  Montgomery  (1995,  p.280)  offer  some  guidelines 
on  the  choice  of  a response  surface  design  (see  also  Box  and  Draper  (1975))-  They 
list  the  following  as  important  characteristics  of  response  surface  designs; 

1.  Result  in  a good  fit  of  the  model  to  the  data 

2.  Give  sufficient  information  to  allow  a teat  for  lack  of  fit. 

3.  -Allow  models  of  increasing  order  to  be  constructed  sequentially. 

4.  Provide  an  estimate  of  'pure'  experimental  error. 

5.  Be  insensitive  (robust)  to  the  presence  of  outliers  in  the  data. 

6.  Be  robust  to  errors  in  control  of  design  levels. 

7.  Be  cost  effective. 

S.  .Allow  for  experiments  to  be  done  in  blocks. 

9.  Provide  a check  on  the  homogeneous  variance  assumption. 

10.  Provide  a good  distribution  of  the  variance  of  tbe  predicted  values. 

There  ate  also  other  criteria  on  which  the  choice  of  designs  is  based.  For  example, 
the  researcher  may  use  the  D,  .4,  E,  C and  J design  optimality  criteria.  However, 
most  of  tbe  results  in  these  design  optimality  criteria  apply  only  to  the  single-response 
case.  .A  few  results  are  also  available  for  muliiresponse  experiments.  Wijesinha  and 
Khuri  (1987a)  presented  results  concerning  multiresponse  D-nptimal  designs  and  Kim 
and  Draper  (1994)  worked  on  the  y-opiimality  criierinn.  .A  description  of  the  results 
obtained  by  Kim  and  Draper  (1994)  and  a generalization  of  their  work  to  other  cases 
are  given  in  Chapter  6. 


1.2  Qbie 


1.  set  up  a scries  of  muUirespoDSe  experiuieots  (design)  that  will  yield  adequate 
and  reliable  predictions  of  the  responses  of  interest  on  the  basis  of  a certain 
hypothesized  modei, 

2.  fit  the  hypothesized  mulllresponse  model  using  data  collected  according  to  the 
design  chosen  in  (1),  and 

that  result  in  the  simuluneous  optimization  of  the  predicted  values  of  the  re- 

The  three  objectives  listed  above  are  addressed  in  Chapters  3 Co  6 of  this  dis- 
sertation. The  problem  of  setting  up  designs,  stated  in  objective  1,  is  tackled  in 
Chapter  6.  A sequential  procedure  for  seeking  a design  that  minimizes  the  integrated 
mean-squared  error  is  presented  in  Chapter  6.  The  second  objective  is  addressed  in 
Chapters  4 and  5.  Tliase  two  chapters  deal  with  the  fitting  of  multiresponse  models 
that  take  into  account  extraneous  sources  of  variation, block  effects.  The  results 
given  are  the  multiresponse  generalizations  of  the  works  done  by  Khuri  (1992, 1994, 
1996b)  In  the  single-response  case.  In  Chapter  3,  a multiresponse  optimization  pro- 
cedure is  introduced  to  obtain  the  simultaneous  optima  of  the  responses  of  interest. 
This  procedure  is  an  extension  of  the  work  done  by  Khuri  and  Conlon  (1981)  to  the 


ISO  where  the 


CHAPTER  2 


REVIEW  OF  LITERATURE 


Response  surface  methodology  (RSM),  as  defined  Myers  and  Montgomery 
(199S,  pT)i  is  a collection  of  statistical  techniques  useful  for  developing,  improv- 
ing, and  optimizing  processes.  These  techniques,  as  described  by  Khuri  and  Cornell 
(1996,  p.3j,  encompass 

I.  selling  up  a series  of  experiments  (designing  a set  of  experiments)  that  will  yield 
adequate  and  reliable  measurementa  of  the  response  of  interest. 

3.  determining  the  optimal  settings  of  the  experimental  factors  that  produce  the 

The  formal  development  of  RSM  started  with  the  work  of  Box  and  Wilson  (1951). 
However,  its  origin  can  be  traced  back  to  the  works  of  .>.  Wishart,  C.  P.  Winsor,  F. 
Yales  and  others  In  the  early  1930s.  RSM  was  initially  developed  for  the  purpcee  of 
determioing  optimum  operating  conditions  In  chemical  industry  applications.  Now.  it 
is  used  in  a variety  of  fields  end  appllcaliuus.  not  only  in  the  physical  and  engineering 


in  the  biological,  clinical, 


Applications  of  RSM  are  mostly  seen  in  situations  where  several  input  variables 
polemially  influence  some  pcrfomtauce  measure  or  characteristic  of  a product  or 
proceaSt  tiiat  is,  input  variables  influencing  one  or  more  response  variables.  Responses 
are  usually  measured  on  a continuous  scale  The  Input  variables,  also  called  control 
or  independent  variables,  are  subject  to  the  control  of  the  experimenter. 

The  basic  response  surface  problem  involves  a response  of  interest,  y' , which  is 

a function  of  a vector  C'  = ((J.Cj' CJ)  variables.  The  relationship  Is 

given  by  S’  = /"(O  +<"  *here  the  form  of  /*  is  unknown  and  <•  is  a term  that 
reprints  other  sources  of  variability  not  accouoted  for  in  /*.  it  is  usually  the  case 
where  is  treated  as  random  error  so  that  it  is  assumed  that  it  follows  the  oormal 
distribution  with  mean  zero  and  variance  cr’.  Thus  E(y')  = r?'  = /'(O-  The  input 

variables  Cj are  said  Co  be  in  their  natural  form  if  tliey  are  expressed  in 

terms  of  the  natural  units  of  measurement,  sucli  as  degrees  Celcius  or  pounds  per 
square  inch,  it  is  convenient  to  analyne  data  without  dealing  with  different  units  of 

measurement  of  the  input  variables.  In  this  case,  Cs Cs  are  transformed  to 

their  coded  form  X|,  xs,  ...,  xs.  These  coded  variables  are  scale  free,  are  centered 
at  Che  origin  of  the  coorilinate  axe.s,  and  have  equal  spreads  in  all  directions  of  the 
coordinate  axes.  In  terms  of  the  coded  variables,  the  true  response  function  can  be 
written  as  >|'  = /'(*)  where  x = (xi.xj xj,). 

Because  the  true  form  of  q*  is  unknown,  it  Is  assumed  that  t)‘  can  be  ^proxi- 
mated  by  a polynomiai  model  In  k variables.  Polynomial  models  derive  their  form  by 
recognising  that  any  continuous,  smooth  functioo  /'  for  which  p = /*(zi,X3,. x^), 
can  be  apyrortmated  locally  about  a point  in  X1.X3, . . . ,Xk  space  with  a Taylor  series 
expansion  of  terms  involving  Xi.xj aud  x^. 

Generally  speaking,  a term  in  the  polynomial  is  said  to  bo  of  degree  d if  it  contains 

the  product  of  d of  the  factors  X|,Z] xt  or  If  d is  equal  Co  the  ' . ,rthe  powers  of 

the  factors  in  the  term.  For  example,  with  three  factors  xi.xo,  and  za,  terms  of  degree 


8 

2arezf,  i|.  i|,  i|ij.  i|ij,  andxjij,  whereas  terna  ofdegreeSarei^,  ij,  ij,  X|ZjXi, 
z}x2.  xfx], . ..,  etc.  polyaomial  is  said  to  be  of  degree  d if  the  model  conlains  terms 
up  to  degree  d.  Hereafter,  when  we  refer  to  a polynomial  model  of  degree  d.  we  shall 
mean  a model  that  containa  all  terms  tn  (he  factors  up  ia  and  including  terms  of 
degree  d.  Polynomial  models  of  degree  1,  2,  and  3 in  the  two  variablea  i|  and  i]  are 
o'  = + diX] -t* 

rf  - /%  + d]i| +d2i}  + 3iiz?  + Aji|  + Aazizj 

rf  = fl>  + dl*l  +A®J  + dl|Z?  + /J!jli  + »JuZ|X2 


respectively,  where  do,  0j,  3j,,  Bjjj,  0,\,  and  djyi,  i.i  = 1,2  ace  unknown  parameters 
or  model  coelhclente  to  be  estimated  from  the  data.  We  also  note  that  a reduced  farm 
of  a polynomial  of  degree  d is  a model  containing  only  some,  or  a subset,  of  tbe  terms 
of  a polynomial  of  degree  d. 

The  response  surface  model  may  be  written  in  tbe  form 

»•  = /■(«, fl*) 


or  equivalently. 


a- = /•(*, 13-) -VC- 


(2-1.1) 


where  /*(x,(3*}  consists  of  all  the  terms  in  the  polynomial  model.  For  a thorough 
discussion  of  R5M,  see  Myers(1971),  Box  and  Draper  (1987),  Khuri  and  Cornell  (1996) 
and  Myers  and  Montgomery  (1995).  \ comprehensive  review  of  this  subject  la  also 
given  by  Myers.  Khun,  and  Carter  (1989). 


2,2.  Multiresponse  Experiments 


In  the  early  development  of  response  surface  methodology,  only  a single  response 
variable  was  considered.  In  many  experimental  situations,  however,  the  experimenter 


» gennally  intetested  io  more  than  one  reapDOse.  For  example,  a chemical  engi- 
neer may  be  icueresied  In  maximizing  the  yield  (response  one)  while  minimizing  the 
cost  [responae  two)  of  a certain  chemical  process.  .-Vlso.  an  agronomist  may  want  to 
maximize  the  harvest  yield  and  minimize  fertilizer  cost  for  a cenain  type  of  crop.  Ex- 
periments in  which  more  than  one  response  is  measured  and  studied  for  a given  setting 
of  the  experimental  factors  (input  variables)  are  called  mu/tiresponse  experiments. 

The  analysis  of  multiresponse  experiments  is  quite  different  from  that  of  single- 
response  experiments.  When  analyzing  data  from  a multiresponse  experiment,  one 
must  take  into  account  possible  interrelationships  that  exist  among  the  responses 
of  interest.  Thus  the  multivariate  nature  of  the  data  obtained  from  multiresponse 
experiments  must  be  recognized  and  considered  in  the  analysis. 

An  important  objective  in  any  multiresponse  experiment  is  the  estimation  of  the 
relationships  between  the  responses  and  the  input  variables,  within  a given  region  of 
interest.  Suppose  that  there  are  r responses  »i,  Sij.  • ••,  Vr.  and  *;  input  variables  x,, 

ij,  ...,  n of  interest.  Let  i = 1.2. u = 1.2 n,  repieseni  n sets  of 

observations  on  each  of  the  r response  variables.  Then  the  relationship  between  the 
responses  and  the  input  variables,  as  given  in  (2.1.1),  can  be  represented  by 

= ■ = 1.2 r;  u = I,2,...,n.  (2.2.1) 

where  Xa  s (x«i.  x^z x^)'  is  the  vector  of  input  variables  at  the  uth  setting,  0, 

is  the  vector  of  unknown  parameters  associated  with  the  ith  response,  <„  is  a random 
error,  and  /*  Is  a respuuse  function  which  can  be  represented  by  a polynomial  of 
degree  d. 

The  estimation  of  the  relationships  in  (2.2.1)  involves  finding  estimates  for  the 
parameters  ft,  i = l,2,...,r.  However,  these  estimates  require  the  collection  of 
observations  for  the  responses  basH  on  given  values  of  the  input  variables.  The  values 
of  Xu.  wherein  the  observed  response  values  will  be  attained,  that  is,  the  design  points, 
are  determined  according  to  an  appropriate  critenon.  Once  a response  surface  design 
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has  been  properly  selected  according  Co  some  criierion,  data  are  collected  and  are 
used  Co  obtain  eslinaces  of  the  unknown  parameters  The  generalized  least  squares 
method  or  the  maximum  likelihood  method  can  be  used  to  obtain  these  estimates  to 
produce  the  fitted  model.  If  the  fitted  models  are  deemed  adequate,  then  they  are 
used  to  determine  the  conditions,  or  settiop  of  the  input  variables,  that  result  in 
estimated  optimal  values  of  the  mean  responses. 


There  are  tvro  approaches  that  have  been  mainly  used  in  the  estimation  of  param- 
eters of  multiresponse  models.  The  first  one  is  the  generalized  least  squares  methud 
introduced  by  Zellner  (1962).  This  approach  uses  his  seemingly  unrelated  regression 
models.  It  also  involves  the  estimation  of  the  variance-covariance  matrix,  £,  asso- 
ciated with  the  r responses.  Zellner  (1962)  showed  Clial  the  parameter  estimates 
obtained  through  this  method  are  at  least  asymptotically  more  efficient  than  those 
obcamed  using  Che  single-equation  application  of  least  squares.  A discussion  of  this 
method  is  given  in  Chapter  3.  The  other  approach,  developed  by  Box  and  Draper 
(1966),  uses  a Bayesian  argument  wherein  non-informative  prior  distributions  are  as- 
signed to  all  the  parameters  iu  the  model,  including  the  variance-covariance  matrix  £. 
This  method  is  briefly  explained  in  Appendix  A.  Note  that  the  Box  and  Draper  (1965) 
approach  is  more  general  than  Zelloer's  (1962)  method  since  it  applies  to  nonflnear 
as  wc)l  as  linear  mulciresponse  models.  Herwevet,  Box  ei  al.  (1973)  showed  that  the 
Former  method  can  lead  to  misleading  results  If  there  are  linear  relationships  among 
the  responses.  Therefore  it  is  important  that  these  dependencies  be  removed  before 
the  Box  and  Draper  (1966)  estimation  procedure  is  used  Box  et  al.  (1973)  proposed 
a procedure  for  detecting  the  existence  of  these  dependencies  when  the  multiresponse 
data  are  subject  to  round-off  errors.  An  extension  of  this  procedure  to  situations  in 


which  the  reepoose  variables  have  different  scales  of  measurement  was  introduced  by 
Khuri  (1990a). 

More  recent  results  on  the  estimation  of  parameters  in  a multiresponse  model  are 
given  by  Kang  and  Bates  (1990).  They  presented  an  approximate  inference  procedure 
which  can  be  used  to  obtain  approximate  cunJideoce  regions  for  the  parameters  in 
the  multlresponse  model.  Their  procedure  mainly  uses  the  estimates  derived  from 
the  Box  and  Draper  (1965)  approach. 

2.4  Multlresponse  Ontimization 

The  literature  addressing  the  problem  of  multiresponse  optimiaation  is  not  as  ex- 
tensive as  that  for  single-response  optimieation.  Khuri  (1906a)  dehnes  multlresponse 
optimization  as  the  determination  of  operating  conditions  on  a set  of  input  variables 
that  yield  optimal  (or  neat  optimal)  values  for  several  response  variables  considered 
simultaneously.  This  particular  problem  of  simultaneous  optimisation  of  responses  is 
addressed  in  the  works  done  by  Harrington  (1965).  Derringer  and  Suich  (1980),  and 
Khuri  and  Coolon  (1981).  See  also,  Khuri  and  Cornell  (1996,  Chapter  7). 

Harrington  (1965),  and  later  Derringer  and  Suich  (1980),  introduced  the  concept  of 
desirability  whereby  each  response  function  is  transformed  to  a desirability  function. 
This  method  make.s  use  of  a deeirability  function  in  which  the  experimenter's  own 
priorities  and  desires  on  the  particular  responses  considered  are  built  into  the  opti- 
mization procedure.  Harrington  (1965)  transformed  each  predicted  response  function, 

to  a desirability  function  d(,  where  0<d(<l.tsl,2,,,,,r,  such  that  d(  in- 
creases as  the  desirability  of  the  conesponding  property  which  y,  measures  increases. 
He  proposed  the  exponential-type  transformation  given  by  d)  = c~i^'',  where  s is  a 
positive  number  specified  by  the  user.  Derringer  and  Suich  (198U)  introduced  more 
geoeral  traDsformatloos  which  offer  the  user  greater  flexibility  io  setting  up  the  de- 
sirability fuDctioDS.  Myers  and  Montgomery  (1995.  pp.  250-251)  gave  the  following 
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description  of  this  approach:  One  and  two-sided  desirability  functions  are  used  de- 
pending on  whether  the  response  is  to  be  minimized  or  maximized  or  has  an  assigned 
target  value.  For  the  ith  response  with  a target  value  level,  .-1.  B,  and  C are  assigned 
so  that  .4  < B < C.  A product  is  considered  unacceptable  if  ^ < 4 or  > C.  The 
value  B is  the  ‘most  desirable  value'  (target  value).  The  quantity  d).  the  desirability, 
is  defined  as 


with  d,*  being  0 if  y,  > C or  p,  < 4.  The  user  then  decides  on  the  values  of  a and 
I according  to  his/her  subjective  impression  of  the  local  desirability  of  Che  response 
under  consideration.  In  the  case  where  the  ith  response  le  to  be  maximized,  the  user 
chooses  a value  B such  that  d|  = 1 for  any  yi  > B.  that  is,  the  researclier  is  satisfied 
if  ^ > B.  However,  the  product  is  assumed  to  be  unacceptable  if  < 4.  In  this  case. 
B = C and  the  desirability  function  is  given  by 


where  4 < y,  < B.  In  the  case  where  the  Ith  response  is  to  be  minimized,  a value  of 
B is  chosen  such  that  ^ < B,  is  considered  desirable  and  produces  a d^  = 1.  A value 
for  M ^ is  considered  unacceptable  and  thus  results  in  a = 0.  Tlie  desirability 
function  is  given  by 


where  B < y,  < C. 

Once  the  desirability  functions  for  all  the  r responses  are  obtained,  a single  com- 
posite response  is  determined.  Tills  single  reeponse  is  the  geometric  mean  of  the  di's 
given  by  D*  s (Ocsi  The  user  then  needs  to  hnd  the  values  of  x iliat  maximize 
O'.  The  optimization  of  O'  is  fairly  simple  since  O*  is  a well-behaved,  continuous 


(2.4.1) 


(2.4,2) 


(2.4.3) 
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functioD.  A value  of  D'  dose  lo  one  Impliea  that  all  responses  are  in  a desirable  range 
sintultaneously.  Note  that  the  geometric  mean  is  only  one  possible  formulation  of  D*. 
Other  functions  of  the  d^’i  may  also  be  used.  Examples  of  these  other  functions  are 
given  by  Myets  and  Montgomery  (199S,  p.253).  For  a more  detailed  discussion  and 
an  example  of  this  approach,  see  Myers  and  Montgomery  (1995.  pp.247-258). 

One  advantage  of  the  Derringer  and  Suich  (19S0)  approach  is  in  the  relative  ease 
In  which  the  optimal  value  of  D'  con  be  obtained.  However,  one  of  its  basic  charac* 
teristics  is  the  subjectivity  on  the  part  of  the  user  in  dehning  the  desirabilities.  This 
is  not  a problem  if  the  user  can  correctly  evaluate  the  desirability  functions.  If  this 
is  not  the  case,  then  it  is  possible  that  the  improperly  ri.sses.sed  desirabilities  will  lead 
to  Inaccurate  results.  Another  disadvantage  of  this  approach  lies  in  the  fact  that  it 
does  not  take  into  account  the  heterogeneity  of  the  variances  of  the  responses  or  any 
correlations  that  might  exist  among  them.  These  are  important  considerations  in  a 
multiresponse  optimization  procedure. 

.Another  multiresponse  optimization  procedure  was  introduced  by  Khuri  and  Con- 
loo  (1981).  This  procedure  is  based  on  the  generalized  distance  approach,  in  this 
procedure,  all  the  response  functions  of  interest  are  assumed  to  depend  on  the  same 
set  of  of  input  variables  and  can  be  represented  by  polynomial  regression  models  of 
the  same  degree  and  form  within  the  experimental  region.  Hence  this  approach  ap- 
plies to  linear  multiresponse  models  of  the  form  y.  = XaJ3,  4-s.,  where  Xo  is  an  n x p 
matrix  of  rank  p and  the  values  of  the  r,'s  have  means  zero,  Vor(e,)  = and 

coo(«„<j)  = (7„I„  i ^ j,  i,j  = 1.2 r.  Let  S be  an  r X r matrix  whose  (t,  j)th 

clement  ie  Og,  t,y  1, 2 r.  Then  an  estimate  of  21  is  given  by 

(i.„i 

where  V = (y,  : v, ; . . . ; y,\.  Without  loss  of  generality  it  can  be  assumed  that  the 
r responses  are  not  linearly  related.  The  prediction  equation  for  the  ith  response  is 
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ji-  = 1 = 1,2 r,  where  *'  (*)  is  a veclor  of  the  same  form  as  a row  of 

the  matrix  Xa  evaluated  at  the  poiot  z,  3,  is  the  OLS  estimator  of  3i,  and  is  equal 

to  i,X',XnY' X'^y,.  1 = 1.2 T (it  can  be  shown  that  this  estimator  coincides 

with  the  generalized  least-squares  estimator  of  A).  1-et  y - - ■yr)'- 

variance-covariance  matrix  of  v(z)  at  a point  (z)  in  the  experimental  region  is  given 
by 

V'nr[i(z)|  = z-'(z)(X'(,.X»)-'z-(z)S.  {2-4.5) 

An  unbiased  estimator,  V'or[y(z)],  of  I'oriviz]]  can  be  obtained  by  using  E in  place 
of  E in  the  above  equation. 

Let  d,  be  the  optimum  value  of  yi  optimized  individually  over  the  experimental  re- 
gion and  let  <d  = (diidr dr)'.  If  these  individual  optima  are  attained  at  the  same 

set,  X.  of  operating  conditions,  then  an  'ideal'  optimum  is  said  to  be  achieved.  How- 
ever. such  an  ideal  optimum  rarely  exists.  In  more  general  situations  one  might  want 
to  consider  finding  compromise  conditions  on  the  input  variables  that  are  somewhat 
favorable  to  all  responses.  This  means  finding  conditions  under  which  the  multire- 
sponse function  deviates  as  little  as  possible  from  the  ideal  optimum,  that  is,  at  such 
conditions  a ‘near’  optimum  is  achieved  for  each  of  the  r predicted  response  functions. 
Such  a deviation  of  the  compromising  conditions  from  the  ideai  optimum  can  be  for- 
mulated by  means  of  a distance  function  which  measures  the  distance  between  ii(z) 
and  Let  this  distance  function  be  represented  by  p|v(z),^).  The  multlresponse 
optimization  then  involves  finding  conditions  on  x that  minimize  this  distance  func- 
tion over  the  experimeotal  region.  Khuri  and  Conlon  (19S1)  gave  several  forms  for 
this  distance  function.  One  of  them  Is  given  by  the  weighted  distance 

P|V(*).«I  = |(v(*)-^)'{I’'ar[y(*)]}’'(jr(z|-  (2.4.6) 

bince  E is  unknown,  V'or(y(z)|  can  be  used  in  place  of  V'ar|v(z)j  in  the  above  equa- 
tion. If  Zq  is  the  point  In  the  experimental  region  at  which  p[y{x),^]  attains  ita 


absolute  miDimum,  and  if  Tn«  is  the  value  of  this  miaiDium,  then  at  zo  the  experi- 
mental conditions  can  be  described  as  being  near  optimal  for  each  of  the  r respoosc 
functions.  The  smaller  the  value  of  mo,  the  closer  these  conditions  are  to  representing 
the  'ideal'  optimum.  A more  detailed  discussion  and  an  example  of  this  approach  is 
given  in  Khuri  and  Cornell  (1996,  pp.966-293). 

An  advantage  oftheKhuri  and  Conlon  (1981)  generalised  distance  procedure  over 
the  Derringer  and  Suich  (1980)  desirability  approach  is  that  the  generalised  dUtance 
approach  takes  into  consideration  the  correlations  among  the  responses  as  wdl  as 
the  magnitudes  of  their  variances.  A disadvantage  of  this  procedure,  however,  is 
the  requirement  that  all  the  response  functions  must  be  represented  by  polynomial 

may  not  be  practical.  This  problem  will  be  addressed  in  Chapter  3. 

Other  optimisation  procedures  that  Involve  several  responses,  but  are  not  con- 
cerned with  'simalUneoue'  optimization,  ere  pven  by  Myers  and  Carter  (1973)  and 
Biles  (1975).  Myers  and  Carter  (1973)  introduixd  the  dual  response  approach  for 
the  optimization  of  one  response  with  constrainUs  placed  on  another  response.  In  the 
dual  response  approach,  the  tvro  responses  are  categorized  as  'primary',  represented 
by  Yf  and  'secondaiy',  represented  by  Y,.  The  goal  of  the  procedure  is  to  determine 
conditions  on  x Chat  optimize  the  primary  response  function  subject  to  au  equality 
coustraint  on  the  secondary  response  function.  The  optimization  problem  is  to 


opSimisB  Y,  subjecl  to  Y,  = c.  (2.4.7) 

order  models  of  the  form 

y,  = b.  + x'b  + x'Bx  + (,  (2.4.8) 

y,  = c.  + x'c  + x'Cx  + i.  (2.4.9) 


where  B apd  C 


k X 1 


5,  and  b„  and 


scalars.  To  restrict  the  search  area  to  a spherical  region  of  radius  p'.  an  additlonai 
constraint  given  by 

x'x=p‘‘  (2.4.10) 

is  considered. 

To  solve  this  optimisation  problem.  Myeta  and  Caner  (1973)  showed  that  a sta- 
tionary point  X has  to  satisfy 

(fl-AiC-AjI)x  = ^i^  (2.4.11) 

where  A|  is  the  Lagrange  multipUet  of  the  constraint  in  (2.4.7),  A:  is  tlie  Lagrange 
multiplier  of  the  constraint  in  (2.4.10).  The  dual  response  approach  consists  of  choos- 
ing Aa  so  that  the  matrix  on  the  left-hand  side  of  the  above  equation  gives  either  a 
maximum  or  a minimum.  This  is  accomplished  by  finding  the  eigenvalues  of  B-A|C. 
The  procedure  calls  for  fixing  the  value  of  A|,  computing  a suitable  value  of  A],  and 
obtaining  X from  (2.4.11).  By  trying  different  values  of  Ai  the  experimenter  obtains  a 
value  of  X that  simultaneously  satisfies  (2.4.10)  and  ensures  a satisfactory  secondary 
response.  Biles  (1973)  extended  these  results  to  the  case  where  there  are  several 
secondary  response  functions. 

Nonlinear  optimiaation  procedures  can  be  used  to  solve  the  dual  response  problem. 
This  is  due  to  the  fact  Chat  usually  both  Che  primary  and  the  secondary  response  func- 
tions are  quadratic  functions.  There  are  several  nunlioear  optimization  procedures 
available  in  the  operations  research  (OR)  literature.  Del  Castillo  and  Montgomery 
(1993)  showed  that  one  particular  nonlinear  programming  algorithm  works  best  in 
solving  the  dual  response  problem.  This  particular  algorithm  is  the  penenlatd  re- 
duced gTTidient  (CRG)  algorithm  . The  following  is  a description  of  the  GRG  method 
(see  Del  Castillo  and  Montgomery  (1993),  p.2Dl):  The  GRG  method  optimizes  a 
nonlinear  objective  function  subject  to  nonlinear  equality  constraints.  Inequality 
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consLraiDts  can  be  transformed  to  equality  coDstraints  ty  introducing  surplus  vari- 
ables. Por  example,  if  a constraint  has  the  form  Xi+X2<  c,  the  surplus  variable  x, 
is  introduced  and  the  constraint  is  rewritten  as  ii  -i-Xn-l-x,  = c.  The  surplus  variable 
is  assigned  a coelEcieal  of  zero  in  the  objective  function,  which  is  the  function  to  be 
optimized. 

Given  a current  feasible  region,  the  GR.G  method  tries  to  find  a direction  of 
movement  that  will  improve  the  value  of  the  objective  function.  Then  it  tries  to  find 
an  appropriate  step  size  for  moving  in  this  direction,  while  checking  to  see  that  the 
constraints  are  not  violated.  The  direction  of  the  movement  is  given  by  the  reduced 
^mdient.  The  negative  of  the  reduced  gradient  gives  the  direction  of  steepest  ascent 
for  a maximization  problem  and  steepest  descent  for  a minimization  problem,  ft  is 
called  a reduced  gradient  because  of  the  way  the  variables  arc  grouped  into  two  sets  by 
tbe  algorithm  at  each  iteration.  One  set  consists  of  basic  (dependent)  variables.  These 
are  the  variables  whose  current  values  are  sinctfy  between  their  bounds.  The  other 
set.  called  the  ooobasic  (independent)  variables,  consists  of  the  remaining  variables 
whose  current  values  are  at  their  bounds. 

If  there  are  m constraints  and  k variables,  usually  m < b,  the  system  of  equations 
generated  by  the  constraints  will  have  more  unknowns  than  equations  and  cannot  be 
solved  explicitly.  However,  the  system  can  be  solved  for  m of  the  basic  variables  in 
terms  of  the  remaining  k-m  nonbasic  variables.  If  the  constraint  functions  arc  'well- 
behaved'.  solving  this  system  can  be  done  by  numerical  methods.  The  solution  of  the 
system  will  be  feasible  and  will  have  at  most  k-m  nonbasic  variables.  The  gradient 
of  the  objective  function  with  respect  tn  the  nonbasic  variables  is  called  the  reduced 
gmdtent  If  the  magnitude  of  the  reduced  gradient  at  the  current  point  is  small,  then 
the  algorithm  stops.  Otherwise,  a line  search  in  the  direction  of  the  negative  of  the 
reduced  gradient  is  performed,  a new  point  is  fonrvl  ^nd  the  process  is  repeated. 
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The  advanisges  of  the  GRG  method  as  discussed  by  Del  Castillo  aod  Montgomery 
(1993)  are  as  follows:  It  allows 

1.  more  coostraints,  heoce  response  functions,  to  be  considered 

2-  the  analyst  to  directly  specify  the  radius  of  the  spherical  region  of  interest  In 
the  dual  response  approach,  rather  than  indirectly  specifying  this  through  the 
Lagrange  multiplier  A|. 

3.  for  other  forms  of  the  constraints  that  define  the  region  of  interest 
FVom  these  advantages  one  can  conclude  that  the  GRG  method  is  more  fiexible  than 
the  Myers  and  Center  (1973)  procedure  and  that  it  allows  for  the  handling  of  more 
then  two  responses.  For  examples  of  the  GRG  approach,  see  Del  Castillo  and  Mont- 
gomery (1993,  pp.  202-204), 

2.5  Blocklna  in  RSM 


Blocking  is  an  experimental  design  technique  which  can  be  used  in  experiments 
that  are  too  large  to  allow  all  runs  to  be  made  under  homogeneous  conditions.  In 
RSM,  there  ore  many  situations  where  this  is  the  cose.  For  example,  a block  may  refer 
to  a group  of  experiments  performed  on  the  same  day  but  several  days  of  experiments 
are  desired,  or  a group  of  experiments  for  which  it  was  possible  to  use  the  same 
batch  of  raw  material  but  several  batches  of  raw  material  are  needed.  Because  of 
the  poesibilitv  of  running  large  experiments  under  non-homogeneous  conditit>ns.  it 
is  important  to  consider  designs  that  can  be  run  in  blacks.  Box  and  Hunter  (19S7) 
indicated  however  that  the  manner  in  which  the  experimental  points  of  a particular 
design  are  assigned  to  the  blocks  greatly  affects  the  efficiency  of  the  estimation  of 
model  parameters  and  the  ease  with  which  these  estimates  are  calculated.  To  solve 
this  problem,  they  introduced  the  concept  of  orthogonal  blocking  for  seennd-order 
designs.  If  a design  blocks  orthogonally,  then  the  inclusion  of  block  effects  in  the 
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model  does  not  affect  the  estimates  of  the  model  parameters  associated  with  the 
model's  polynomial  effects.  The  intercept  estimate,  however,  depends  on  the  blocltB. 
The  analysis  of  the  experiment  then  proceeds  exactly  as  if  there  were  no  biocking 
except  for  a reduction  in  the  residual  snm  of  squares  due  to  the  removai  of  the  blocks 
sums  of  squares.  Box  and  Hunter  (1957,  1961a.b)  developed  the  general  conditions 
that  give  rise  to  orthogonal  blocking  in  the  second-order  case.  They  also  presented 
orthogonal  blocking  arrangements  for  some  central  composite  designs  (CCD),  Box 
and  Hunter  (1957)  noted  that  one  aspect  which  mekes  these  orthogonal  blocking 
arrangements  particularly  useful  arises  out  of  the  nature  in  which  CCDs  are  used,  that 
is,  in  the  exploratinn  of  response  surfaces,  trials  are  usually  performed  sequentially 
and  often  have  as  their  objective  the  optimization  of  some  response.  For  example, 
suppose  that  four  variables  were  studied  with  the  intention  of  finding  an  optimum 
set  of  conditions  for  a new  process  using  the  method  of  steepest  ascent.  The  method 
of  steepest  ascent  is  an  optimization  technique  which  involves  sequential  movement 
from  one  region  to  another,  to  search  for  conditions  in  which  the  process  output  value 
is  improved.  In  this  case,  a half  replicate  of  a 2*  factorial  design  with  two  center 
point  replicates  (block  1)  could  first  be  run.  This  design  vrould  supply  estimates  uf 
the  four  hrst-order  terms,  combinations  of  the  interaction  terms  and  of  the  sum  irf 
the  four  quadratic  terms.  If  the  values  of  the  coefficient  estimates  of  the  htst-order 
terms  were  very  large  compared  to  their  siaodacd  errors,  progress  could  probably 
be  made  without  a more  elaborate  design  by  moving  in  the  indicated  direction  of 
iocreasing  yield.  Continuous  movement  would  be  made  until  improvement  in  that 
direction  was  exhausted,  .kt  this  point  a new  ftrst-order  design  would  be  started  at 
the  improved  set  of  conditions.  If  the  first  order  terms  were  small,  ora  tentative  move 
bad  not  met  with  success,  the  second  half  of  the  factorial  design  together  with  two 
additional  center  point  replicates  (block  2}  could  be  carried  out,  and  the  situatiun 
is  reviewed  In  the  light  uf  the  set  of  20  trials  which  are  now  available.  Finally,  if 


evfdeace  indicstes  that  further  progress  can  be  attained  only  by  fitting  the  second* 
degree  equation,  the  third  bfock  of  ten  points  consisting  of  the  dght  axial  points  and 
two  center  points  would  be  added.  The  complete  set  of  thirty  trials  could  theo  be 
used  to  provide  efficient  estimates  of  the  beet  fitting  second-order  model  and  further 
progress  would  follow.  Thus,  bistorically,  blocfclng  in  RSM  is  mainly  used  as  a tool 
in  sequential  experimentation,  fbr  a more  detailed  discussion  of  orthogonal  blocking 
for  second-order  designs,  see  Box  and  Draper  (1987,  pp.509-520),  Khuri  and  Cornell 
(1998,  Chapter  8)  and  Myers  and  Montgomery  (1998,  pp.  328-341).  Results  on  non- 
orthngonal  blocking  for  second-order  models  are  also  gsven  by  Dey  and  Das  (1970) 
and  Khuri  and  Cornell  (1998,  Chapter  8). 

Some  of  the  more  current  results  on  blocking  in  RSM  arc  given  by  Cook  and 
Nactsheim  (1989),  Adhikaiy  and  Panda  (1990),  Khuri  (1992, 1994, 1996b)  and  Khuri 
and  Cornell  (1996,  Chapter  8).  Cook  and  Nacisheim  (1989)  introduced  computer- 
aided  blocking  of  factorial  and  response-surface  designs  while  Adhikaty  and  Panda 
(1990)  presented  a sequential  method  fur  constructing  second-order  rotatable  designs 
in  uon-orthogonal  blocks.  Khuri(L992)  considered  reeponse-surface  models  with  ran- 
dom block  effects  and  generalised  the  result  given  by  Box  and  Hunter  (1987)  ou 
the  general  conditions  for  response  surface  designa  to  block  orthogonally.  Details  of 
his  work  arc  given  in  Chapter  S.  Khuri  (1994)  also  presented  results  on  the  effect 
of  blocking  on  the  estimation  of  the  mean  response,  on  the  prediction  variance  and 
on  the  optimum  of  the  response  surface.  These  results  are  described  in  Chapter  4. 
More  recently,  Kliuri  (1996b)  described  results  on  response  surface  models  with  mixed 
effects.  A descriptiou  of  these  are  also  given  in  Chapter  5. 

2.6  Multiresnonse  Design  Criteria 

In  selecting  a design,  most  researchers  rely  on  a design  optimality  criterion.  Mrjst 
design  optimality  criteria  for  single-response  experiments  ace  bssed  on  some 
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properties  of  the  X X matrix  for  a design  matrix  X One  such  rriuirion  is  the  D- 
optimality  criterion.  D-npiimeliiy  is  one  of  llie  most  pt^ular  design  criterion.  D- 
optimal  designs  minimise  the  determinant  of  (X'X)''.  Its  appeal  is  based  on  the 
property  that  the  determinant  of  X'X  is  inversely  proportional  to  the  square  of  the 
volume  of  the  confidence  region  (hyperellipsoid)  on  the  model  parameters.  Thus  a 
large  determinant  of  X'X,  or  equivalently  a small  determinant  of  (X'X)  implies 
better  estimation  of  the  model  parameters. 

In  the  multiresponse  case,  one  needs  to  consider  criteria  that  involve  all  the  re- 
sponse variables.  Onesuchcriterion  was  developed  by  Fedorov  (1972).  Fedorov  (1972) 
introduced  a sequential  procedure  that  cau  be  used  to  obtain  multiresponse  D-optimal 
designs.  However,  his  procedure  requires  knowing  E.  the  variance-covariance  matrix 
associated  with  the  r responses.  Wijesinha  and  Khuri  (19S7a)  modified  Fedorov's 
procedure  by  replacing  £ with  Zellner's  (19B2)  estimate.  .Another  design  criterion 
presented  by  Wijesinha  and  Khuri  (1987b),  is  one  based  on  maximising  the  power 
of  the  multivariate  lack  of  fit  test  introduced  by  Khuri  (1986).  This  criterion  is  an 
extension  of  the  method  used  by  Jones  and  Mitchell  (1978)  in  the  single-response 
case.  Khuri  (1990b)  introduced  the  concept  of  multiresponse  rotatability,  based  on 
the  work  done  by  Box  and  Hunter  (1957).  The  concept  of  rotatability  was  first  intro- 
duced by  Box  and  Hunter  (1957)  and  has  since  become  one  of  the  more  popular  and 
important  design  criteria.  With  a rotatable  design,  the  prediction  variance  at  a point 
X is  a function  only  of  the  distance  between  the  point  x and  the  center  of  the  design. 
Tills  implies  that  with  a rotatable  design,  the  prediction  variance  is  the  same  at  all 
points  X that  are  equidistant  from  the  design  center.  Khuri  (1990b)  extended  this 
property  to  the  multiresponse  case.  He  defined  multiresponse  rotatability  as  a design 
property  under  which  the  variances  and  covariances  among  the  predicted  responses 
are  constant  at  all  points  that  are  equidistant  from  the  design  center.  H*-  showed  that 
multiresponse  routabllity  can  be  achieved  if  and  only  if  the  design  is  rotatable  for 


a single-respoDse  model  that,  haa  Ihe  highest  degree  among  all  the  response  models 
under  coosideration. 

In  choosing  a response  surface  design,  an  cxperiuientcr  may  want  to  consider  a 
design  criterion  that  takes  into  account  both  the  bias  and  the  variance  of  the  fitted 
model.  The  integrated  mean-squared  error  (IMSE)  criterion,  introduced  by  Bo*  and 
Draper  (1959),  considers  both  properties.  The  IMSE  criterion  is  also  known  as  the 
J-criterion. 

The  J-  criterion  is  described  as  follows  (see  Khuri  and  Cornell  (1996),  Chapter 
6)r  Suppose  that  the  fitted  model  for  the  mean  response  is  given  by  II  = + « 

of  order  d|.  Let  the  true  model  for  the  mean  response  be?)  = X\0i  + X‘,0j  of  order 
di  (dj  > d| ).  Then,  at  each  point  i in  the  experimental  region  »,  the  mean  squared 
error  (MSE)  of  j(*)  is  .V/S£[p(r))  = £[!i(j:)->/(x)J’  where  j(*)  is  the  least  squares 
estimate  of  the  mean  response  value,  h(x),  at  x obtained  with  the  fitted  model.  If 
this  quantity  is  averaged  over  the  region  S,  we  obtain;  n/,£|v(x)-  >;(i)|*di  where 
f!*'  = /qdx.  Novi,  in  order  to  be  able  to  center  our  attention  on  designs  which  may 
or  may  not  contain  a fixed  number  of  points  and  also  to  remove  the  dependence  of 
the  mean  squared  error  of  p(x)  on  the  value  of  the  variance  the  criterion  is  put 
on  a ‘per  observation'  basis  and  is  expressed  as  d = {Sil/a^)St>E{v{x)  - i)(x)|“dx. 
The  difference  j(x)  - q(x)  can  be  partitioned  as;  y(x)  - t|(x)  = {y(x)  - £(v(x)]}  + 
(£[i)(xj|  - h(x)}.  This  partitioning  enables  us  to  separate  J into  ivro  parts,  namely, 

J = (.Vn/er’)^KQr(s(x)]dx-h(.Vn/cr’)jr{£[y(x)| -q(x))’dx 
= V + B 

The  quantities  V and  B are  called  the  average  variance  of  y(x)  and  the  average 
squared  bias  of  y(x),  respectively  The  ./-criterion  is  concerned  with  minimiring  V 
and  8.  but  regards  the  minimisation  of  B to  be  more  Important.  Kiefer  (1972) 
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introduced  a multiresponse  version  of  the  integrated  mean-squared  error  (IMSE)  cri- 
terion. He  presented  an  iterative  procedure  for  minimizing  the  LMSE  based  on  the 
results  given  by  Karsoo,  Manson  and  Hader  (1969),  under  the  assumption  that  the 
variance-covariance  matrix,  E.  associated  with  the  r responses  is  known.  More  re- 
cently, Kim  and  Draper  (1994)  presented  results  for  minimizing  the  IMSE  for  two 
responses,  assuming  aiso  that  the  variance-covariance  matrix  E is  known.  A more 
detailed  discussion  of  this  criterion  is  given  in  Chapter  6.  Examples  and  applications 
of  the  y-criterion  can  be  found  in  Khuri  and  Cornell  (1996,  pp.207-231)  and  Myers 
and  Montgomery  (1995,  pp.406-421). 


CHAPTER  3 


OPTIMIZATION  OF  MULTIRESPONSE  FUNCTIONS 


Oae  of  the  objectives  of  response  surface  methodology  is  the  detecminatioo  of 
conditions  on  the  input  variables  that  produce  optimum  values  of  the  predicted  re- 
sponse. In  this  chapter  we  present  a method  that  can  be  used  to  find  these  optimum 

mulliresponse  optimisation  can  be  seen  in  many  practical  situations.  For  example, 
in  the  pharmaceutical  area,  the  sucntist  is  not  only  concerned  with  maximisation 
of  the  efficacy  of  a drug  but  also  with  the  minimlaatlon  of  undesirable  side  effects 
of  the  drug.  The  optimisation  of  a multirespouse  function  must  take  into  account 
any  interrelationships  that  exist  among  the  responses  of  interest.  Failure  to  do  so, 
may  lead  to  optimal  conditions  for  one  response  that  are  not  optimal  for  the  other 
responses. 

The  simultaneous  optimisation  of  several  response  functions  was  addressed  by 
Khuri  and  Conlon  (1981)  using  what  is  called  the  generalised  distance  approach.  A 
description  of  this  approach  was  given  in  Chapter  2.  In  this  opiimisaiion  procedure. 
Khuri  and  Conion  (1981)  considered  the  case  in  which  the  responses  are  represented 

case  that  models  of  the  same  degree  and  form  arc  appropriate  for  all  the  responses. 
In  some  instances,  the  design  chosen  to  run  the  multiresponse  experiment  may  not 


provide  s sufficient  number  of  runs  or  points  to  support  fitting  a higher  degree  model 
that  some  of  tiie  responses  {but  not  all)  demand.  U is  also  possible  that  the  'best' 
model  for  one  response  may  not  be  the  'best'  for  the  other  responses.  Thus  we  extend 
the  generalized  distance  approach  presented  by  Khuri  and  Conlon  11981)  to  the  case 
where  the  polynomial  regression  models  are  not  necsssarily  of  the  same  degree  or 


In  this  section  we  will  establish  a single  model,  that  is,  a multiresponse  model, 
which  will  represent  all  the  responses  of  interest.  The  true  form  of  each  response 
function  Is  assumed  to  be  adequately  approximated  by  a polynomial  model.  Once 
the  multiresponse  model  is  established,  estimates  of  the  model  parameters  are  derived 
by  using  the  generalized  least-squares  method. 

Suppose  that  there  are  r responses  of  interest,  namely,  Vi.  Pi.-  . . Vr>  each  of  which 
depends  on  the  same  set  of  1:  independent  variables  denoted  by  Xi.Z;, ...  within 
an  experimental  region  Ti.  Using  the  results  given  by  Khuii  [1990b),  the  model  for 
the  ith  response  is 


where  Pt(x)  is  a polynomial  of  degree  d,  defined  at  x = . , , , Ze)'  and  Ci  is  a 

random  error.  Model  (3.2.1)  can  be  expressed  as 


where  /[(x)  is  a vector  of  ^ elements  consisting  of  powers  and  cross-products  of 

powers  of  Z|,Z] Zk  up  to  degree  d,  and  whose  first  idemeot  is  equal  to  one.  and 

is  a vector  of  unknown  parameters  of  order  Pt  x 1. 


«.(x)  = pi(x)  + <i. 


i=  1,2, 


(3.Z1) 


».(»)  = /:(»)7.-br.,  1 = 1.2 r. 


(3.2.2) 
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If  observations  are  made  on  the  r responses  at  n design  settings,  then  from  (3.2.2) 
we  obtain  the  modeb 

V,  = t/,7,  + fe,  1 = 1,2 r (3.2.3) 

where  y,  and  e,  ate  the  vectors  of  observations  and  random  errors,  respectively,  for 
the  ith  response,  and  U,  is  of  order  n x e,  and  rank  whose  first  column  is  1„.  The 
modeb  in  (3.2.3)  can  he  expressed  as  a single  multiresponse  model  of  the  form 

V = C/7  + e (3.2.4) 

where  y = (u'l.pi y',)'.  7 = (Vi.Ti.  - ,7r)’'  ^ = (iiag{U,.U-i U,). 

The  variance-covariance  matrix  of  e is  0 In,  where  £ = (oty)  b the  variance- 
covariance  matrix  of  the  r responses.  The  best  linear  unbiased  estimator  of  7 is 

■y  = (ir(E-‘  0 /„)t/]-'l/'(E-'  0 /„)y.  (3.2.5) 

The  variance-covariance  matrix  of  -y  b 

t'or(7)  = (U’(E-’  0 In)U]-'.  (3.2.6) 

Let  y[x)  be  the  rxl  vector  of  predicted  responses  at  a point  z in  the  experimental 
region;  that  is, 

jr(z)  = [iii(z),j3(*) »-(»)]'•  (3.2.7) 

where  yi(z)  = /|(z)7,  and  7,  b the  best  linear  unbiased  estimator  of  7,  obtained 
from  (3.2.5).  Formula  (3.2.7)  can  then  be  written  os 

y(i)  = A'(x)7,  (3,2.8) 

where  A'(z)  = dias(/'|(i), /',(z) /i(i)),and  7 = (7’n -yi. . . - ,7^)'.  From  (3.2,6) 

and  (3.2.8)  the  variance-covariance  matrix  o(  y[z)  is  given 


VV(ir(z)|  = A'(z)|l/’(E-'  0/„)[f]-'A(i). 


(3.2.9) 
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The  matrix  S a usually  unknown.  Tocstimarethe  elements  of  S we  use  the  consisteni 
estimator  given  by  Zellner  (1962)  (see  also  Srivastava  and  Giles  (1987,  pp.  13-17)]. 
This  estimator  is  S = (sy)  where 

Si,  = (3,2.10) 

Note  that  other  estimators  of  S may  be  used.  A brief  description  of  the  maximum 
likelihood  rstimator  of  £ is  given  in  .Appendix  .A.  However,  this  approadi  was  devel- 
oped under  the  Bayesian  framework  and  also  under  the  assumption  that  e fdlows  a 
multivariate  normal  distribution.  The  assumption  of  normality  for  e in  our  case  is 
extraneous  in  the  sense  that  we  do  not  need  it  for  the  optimization  procedure  that 
will  be  presented  in  the  next  section.  Hence,  for  our  optimization  procedure,  wo  will 
use  S to  estimate  £. 

Using  5 given  In  (3.2.10).  7 is  estimated 

7,  = [U'(S-'  a/„)y  (3,2.11) 

and  the  corresponding  vector  of  predicted  responses  is  ^veo  by 

y,(x)  = A'(x)7,  (3.2.12) 

Thus,  we  approximately  have 

for(»,(i))  = A'(l)(l/'(S-'  »I„)t/|-'A(x).  (3.2.13) 

The  assumed  uonsingularity  of  £ does  not  ensure  the  nonslngularity  of  5.  Srivastava 
and  Giles  (1987.  p.l6)  showed  that  if  the  number  of  responses,  r,  is  greater  than  the 
number  of  observations,  n,  then  S will  be  singular  with  probability  one.  They  also 
showfd  that  r < n is  a necessary  but  not  sufRcieut  condition  for  the  nonsingularity 


of  S. 


3.3  Proetdufs  for  ihf  Simuitt 


Qptimi3 


nf  1>1>  RMiyinoti 


Id  ihc  previous  section  we  have  obtained  the  estioialors  of  the  model  parameieni. 
This  allows  us  to  proceed  with  the  simultaneous  oplimiaotion  of  the  responses  via  the 
geoeralized  distance  approach. 

Multiresponse  upiimiratioo  under  the  generalized  distance  approacli  amuums  to 
finding  conditions  on  x that  minimize  the  distance  between  the  predicted  responses 
and  their  corresponding  optimum  values.  Let  be  the  optimum  value  of  j,,(x) 
optimized  individuallv  over  an  experimental  region.  where  is  the  ith  element 
of  yj{*)  in  (3.2.12).  Lei  = (di,  ©e, . . . , Also,  lei  p1v,(»),  be  a distance 
function  which  measures  the  distance  of  y,[x)  from  <t-  Consider  the  distance  function 
given  by 

which  can  also  be  writteo  as 

p[v.(»).0l  = ((ir,(x)-#)'lA'(x)(Cr(S-'3X„)Cf)-'A(a:)i-‘(ir,(ic)-d))}''^  (3.3.1) 

To  minimize  p[y,|z),0),  we  need  to  take  the  variability  of  which  is  a random 
vector,  into  account.  Tb  do  so,  we  apply  the  following  procedure  given  by  Khuri  and 
Conlon  (1981):  Let  C be  the  optimum  value  of  the  true  mean  of  the  ith  response. 

optimized  iodividually  over  Tt  and  let  ( = (Cii^ Cr]''  objective  la  to  find 

ao  z in  IZ  such  that  p|v>(x),C|  has  an  absolute  minimum  over  Ti.  However,  p is  a 
function  of  ( which  is  unknown.  As  a result,  it  is  impossible  to  mioimize  p directly. 
In  this  case,  we  will  minimise  an  upper  bound  instead.  This  upper  bound  can  be 
established  by  setting  up  a confidence  region  on  C-  bet  be  a confidence  region  for 
C,  then  whenever  C € 

p[Ve(»).Cl  < '"‘“ipei'^bii'rCxI.vi.  (3.3.2) 


Since  C i*  unknown,  vre  adopt  a conservative  approach  and  minimise  the  function  on 
the  right-hand  side  of  the  above  inequality  over  The  following  result  is  needed  for 
this  approach: 


Thmrm  ^..1.1  The  confidence  regions  sioen  by 

9i,<<i<9!.,  i = 1.2 r (3.3.3) 

hold  simultoneously  to  produce  a rectangular  confidence  region  with  an  approximate 
ron^cnce  coc^ctcnr  of  at  least  I - o'.  = I - ra.  nJienr 

SI.  = A - C,lU,LS)c„ii.  92.  = ®.-^  (3.3.4) 

In  (S.S.i).  C^lU.t.S)  = l/;(€.)V"/.(i)l‘'’.  ai/iere  {.  u the  pomi  si  loAicIi  yu(x) 
sttoina  Its  mdividval  optimum  *.  /[({,)  « the  oeetor  f\  etisiuo(erf  si  j,.  v“  u the 
ith  diagonal  .ruhmstnz  of 

(s"U\V,  3'tU',Ut  ...  s'-t^itf, 
s”U,U2  ...  s^U',U. 

s’'U‘,U,  s'‘U',l/.j  ...  s"U',U, 

aads"  IS  the  ijth  element  of  S~'.  This  rectangular  confidence  region  can  be  used  in 
place  ofV(  in  fS.S.Sf. 

Proof:  Let  denote  the  value  of  x at  which  the  true  mean  of  the  ith  response 
achieves  an  optimum  C aver  K.  ignoring  the  remaining  responses.  An  unbiased  esti- 
mator of  C.  is 

C = /;«.)7„  ■ = 1,2 r. 

where  /(({,)  is  the  vector  /( , valu  .ted  at  4,,  for  which  the  variance  is 


V'ar(C)  = /:(«,)V" /,(«,),  . = I, 


where  V”  ia  the  ilh  diagonal  submatrix  of  |CT[E  ^ written  in  the  form 


a"U\Ux 

a^'V'^V, 


a^’V\V. 

n^'U\V, 


V n'‘U',V, 


r.  = r. 

I/:(€.)V  /,(«.)!* 

where  v“  is  iheiih  diagonal  submatrix  of  [C/'(S’' (t,  = /{(€,)7i,.  and  7^ 
is  obtained  from  equation  (3.2.11).  Then,  using  the  results  of  Zellner  (1962,  p.353), 
7,  has  the  standard  normal  or  Z-distrihution  asymptotically.  .'Vow,  P(ird  < i,/j)  is 
approximately  equal  to  1 -a,  i = 1, 2. . . . , r.  Using  the  Bonferroni  inequality,  we  have 
approximately 

If  we  let  1 - ra  = 1 - o',  then  the  inequalities 

w„<C,<iio„  1 = 1,2 r 

hold  simultaneously  with  an  approximate  confidence  eoelBcieni  of  at  least  1 — u*. 
where 

= i.  - i/;(«,)V''7,(«,)i''’=,/s. 
sq.  = Cw  + l/;«.)V'7,(«.)l‘'’--s/r. 

If  is  replaced  by  then  the  above  inequalities  become  inequalities  (3.3.3).  It 

follows  that  an  approximate  confidence  region  for  C = ((i,C2 Cr)‘  with  confidence 

coefficient  of  at  least  1 — a*  is  given  by  the  rectangular  region  determined  by  inequal- 
ities (3.3.3). 


1b  this  sectioD,  two  examples,  each  illusttating  the  direct  minimization 
distance  function  p,  are  discussed.  The  first  example  uses  the  results  by  Tseo  et  al, 
(1983)  (see  also  Khnri  and  Cornell,  1996,  p.295)  and  deals  with  four  response  variables 
and  three  input  variables.  The  second  example  uses  the  study  made  by  [lichen  ei  al. 
(1974)  which  consists  of  three  response  variables  and  five  input  variables.  For  these 
two  examples  we  show  that  the  'best'  models  in  terms  of  the  adjusted  R^,  and 
the  Cp  criteria  are  not  the  same  for  all  the  responses  Hence  the  procedure  described 
In  Section  3.3  may  be  applied  to  determine  the  simultaneous  optima  of  the  responses. 

Exampie  1:  Tseo  et  at.  (1983)  investigated  the  effects  of  washing  temperature 
(ii),  washing  time  (xj),  and  washing  ratio  (ij)  on  springiness  (pi).  thiobsrbituric  acid 
number  (w).  percent  cooking  loss  (yj),  and  whiteness  index  (jm)  for  minced  mullet 
flesh.  It  is  of  Interest  to  simultaneously  maximise  pi  and  y,  and  minimize  yy  and  ps. 
The  design  settings  in  the  original  and  coded  variables  and  the  multiresponse  data 
ate  given  in  Table  3.1. 

Before  we  can  find  the  simultaneous  optima  of  the  four  respooses,  we  Deed  to  find 
the  ‘best’  multiresponse  model  for  the  data.  To  do  so,  we  can  use  the  PRESS  aod  the 
Cp  statistics  for  seemingly  unrelated  regression  models  derived  by  Sparks  (1967).  As 
an  additional  tool  <we  can  also  use  a measure  of  goodness  of  fit  for  seemingly  unrelated 
regressions  given  by  .VlcElrc^  (1977).  namely.  McElroy's  R^.  Detailed  discussions  of 
the  PRESS  and  Cp  criteria  for  multiresponse  models  and  McElroy's  R^  are  given 
in  Appendix  B,  Seciitms  B.l  and  B.2.  respectively.  Because  of  the  large  number  of 
subset  models  that  we  need  to  consider  here,  we  will  follow  the  recommendation  of 
Sparks  (1987)  for  reducing  the  number  of  subset  models  considered:  First,  by  using 
I'le  R^,  adjusted  R,  and  the  Cp  criteria,  we  find  all  the  ‘good’  subset  models  for 
for  t = 1, 2,3,4.  For  each  pi,  i = 1, 2.3,4,  we  start  with  the  model  where  al)  the  linear 
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Table  3.1.  Experimemal  Deaipi  aud  Respopse  Values  (exampWj__ 

Original  Input  Variables  Coded  Input  Vsnabla Retponses 

29.31  29.30  50.36 


1.982  1.90  20-16 


Source:  Khun  and  Cornell  (1996.  pp.296-296) 


' 30.19 
I 50.86 
I 50.84 
I 30.93 


terms  ii.ir,  and  ij  are  included,  then  we  increase  the  number  of  terms  in  the  model 
by  adding  the  second-degree  (puce  and  cross-product)  terms  one.  two.  or  three  at  a 
time  and  then  stop  when  all  the  terms  in  the  model  form  the  complete  second-degree 
model.  Under  the  Cp  criterion,  a subset  model  is  considered  'good'  if  its  Cp-value 
is  dose  to  p,  the  number  of  paiameteis  in  the  subset  model  being  considered.  Using 
R‘  and  adjusted  i?,  we  say  that  a subset  model  is  'good'  if  the  iP  and  adjusted  fl* 
\-alues  are  close  to  1,  We  consider  all  of  tliese  criteria  to  determine  the  'best'  models 
for  each  of  the  four  responses  of  interest.  .411  these  quantities  were  obtained  by  using 
the  procedure  PROC  REG  of  SAS.  The  models  that  are  chosen  for  the  four  responses 
along  with  their  A’,  adjusted  and  Cp  values  are  given  in  Tables  3.2  - 3.5. 

Second,  we  examine  all  possible  multiresponse  models  of  the  form  ^ven  in  equa- 
tion (3.2.4)  by  includin'  u its  ith  equation  a good'  subset  model  for  y,,  i = 1. 2, 3, 4. 
This  means  that  w«  consider  all  possible  combinations  of  all  the  models  for  yi,}n,!a. 


34 


and  V,  given  in  Tables  3.2  - 3.5.  For  example,  uang  Tables  3.2  • 3.5,  one  possi- 
ble muliiresponse  model  is  formed  by  using  model  (3.4.1)  for  i/i,  model  (3-4-5)  for 
n,  mode!  (3-4.8)  for  jn,  and  model  (3.4.10)  for  y,.  The  next  passible  multiresponse 
model  is  given  by  the  coniblnation  of  models  (3-4-1).  (3-4.5),  (3.4.8),  and  (3.4.11), 
then  the  next  set  is  given  by  models  (3.4.1),  (3.4.5),  (3-4.8),  and  (3.4.12),  and  so  on 
until  we  get  the  last  set  given  by  models  (3.4.4),  (3.4.7),  (3.4.9),  and  (3.4.12).  If  we 
do  this  vee  should  get  a total  of4x3x2x3  = 72  possible  multiresponse  models.  For 
each  resulting  mnitiresponse  model.  Sparks'  PRESS  statistic  and  McElroy's  B}  are 
then  computed.  The  results  for  the  72  multiresponse  models  are  given  in  Table  3 6. 
The  Cp  suti-uic  was  not  computed  for  the  resulting  multiresponse  models  because 
according  to  Sparks  (1987),  the  Cp  criterion  is  only  useful  if  (n  - i:,’'.,  ?,)  is  fairly 
large,  that  is.  greater  than  25.  In  this  case,  the  largest  value  that  this  quantity  will 
iiave  is  17  - 8 = 9,  which  is  quite  small. 

Finally,  from  Thble  3.6  we  choose  the  'best’  multiresponse  model,  that  is,  the 
model  with  the  smallest  PRESS  statistic  value  and  R’  value  close  to  1.  Wo  select 
model  number  48  in  Table  3.6  because  it  clearly  has  the  smallest  PRESS  sutistic 
value  and  it  also  has  an  value  quite  close  to  that  of  the  highest  value  associated 
with  model  number  1 in  Table  3.6.  We  can  then  use  model  (3.4.4)  for  y\,  model 
(3.4.7)  for  <n.  model  (3.4.9)  for  ys,  aod  model  (3.4.11)  for  y,. 

Next,  we  need  to  determine  the  generalized  least  squares  estimate  7 for  the  mul- 
tiresponse model  formed  by  using  models  (3.4.4),  (3.4.7),  (3-4-9),  and  (3.4.11).  The 
S matrix  is  given  in  Table  3.7.  Using  thsse  values  the  matrix  7,  is  computed  and  the 
resulting  values  for  7ier72e,  7je,  and  7^,  are  as  follows: 

fi,  = (1.8766,-0.0974,-0.0009.0.0091,-0.1020,0.0044,0.0037)', 

7i,  = (22.6128,5.6609,-0.1719,-1.2268,7.8551,0.0663,2.6920.0.1752)'. 

7^,  = (17.9782,0.7442.-0.0120,-1.0710,3.4414,0.7004,1.6328,1.3020,1.9716)', 


72  3:4.!  tl?  IaI  3I12  ameisie  S 
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r.hl,3.7.  Th^Smtlrli 
0.001288  0.001S32  -0.024424  -0.067885 

0.001288  0.001632  -0.024424  -0.067885 

-0-024424  -0.397570  1.990787  2.789313 

■0.067885  1.978386  2.789313  11.090586 

= (51.9454,2.4364.0.6307,0.2692,-3.43691.-0.0358,-1.2053,0.6016,0.1260)'. 
Thu6  the  models  for  the  four  responses  are  given  by 

yi,  = 1.8766  - 0.0974H -0.0009I^-^-00091Z5  -0.1020I|-l- 
0.0044I|  -h  0-0037i| 

yu  = 22.6128  4-5-6609i,-0.1719rj-  1.2268*5  4- 7.8551iJ-f 
0.0663*5  + 2.6920iin  4-  D.i752ij*5 
W,  = 17.9782-t-0.7442i| -0.0120*5-  10710*s-l-3.4414i;-f 
0.7004*1  -H  1.6328*5  + 1-3020*1*5  + 1.9716*1*3 
ji,  = 51.9454  + 2.4364*i-e0.6307i2-hO-2692*3  - 3.43691*?- 

0.0358*5  - 1.2953*1*5  -I- 0.6010*1*3 -I- 0.1260*3*3  (3-4.13) 

Using  the  y,,  values,  we  can  now  find  the  individual  and  simultaneous  optima 
uf  the  four  respooses.  The  individual  optima  were  obtained  by  using  a FORTRAN 
program  written  by  Conlon  (1992)  which  applies  the  optimization  algorithm  of  Price 
(1977).  .A  description  of  this  algorithm  is  given  in  Appendix  B.  Section  3.3.  The 
individual  optima  and  their  locations  are  given  in  Table  3.8.  .Note  that  the  experi- 
mental region  is  given  by  Z,’=i  i?  < 3.  The  distance  measure  o in  equstion  (3J.1) 
is  directly  minimized  using  the  procedure  PROC  IML  of  SAS  where  16,000  values 
of  p were  computed.  The  minimum  value  of  p,  the  simultaneous  optima  of  the  four 

tuna  of  yii^iVs,  and  are  close  to  their  corresponding  individual  maxima/minlma. 


T»ble  3.8.  The  Individutl  OptimB  of  the  Four  Responses 

ReapoDM  Individml  OpLlma Location 

in,  01=1.9271  (-O.M28. -0.0113, 1.66901' 

jj,  03  = 18.8366  (-0.5524,1.3082,0.9917)' 

jv  01  = 17.5833  (-0.3096,0.3028.0-5161)' 

0<  = 53.4116  (0.5488, 0.S427, 1.5503)' 


Table  3.9^h^imultan»u^gtta^^h^bu^eaponse9 
ResooDM  Sitnullanecnu  Optima 

Vi,  19166 

V3,  18.9631 

in,  18.4019 

V„ 51.9021 

Location  x = (-0.4917,1.4904,0.7079)' 
min  p 1.8514 


U should  also  be  noted  that  the  minimax  approach  based  on  Theorem  3.3.1  can  also 
be  applied  in  this  example. 

The  optimiaation  procedures  using  the  desirability  function  aproach  by  Derringer 
and  Suich  (1980)  and  the  nonlinear  optimizatioo  procedure  using  the  GRG  algorithm 
were  also  used  in  this  example.  These  two  procedures  were  discussed  in  Chapter  2. 
Section  2.4.  The  results  for  the  desirability  function  approach  are  given  in  Tbble  3.10 
and  the  results  for  the  GRG  approach  are  given  in  Table  3.11.  The  simultaneous 
optima  and  the  desirability  values  in  Table  3.10  were  computed  using  the  Design- 
Expert  software  (see  Design-Expert  Version  4.0  User's  Guide  for  details  on  how  to 
use  the  program)  and  the  computer  output  is  given  in  Appendix  C.  When  using  the 
desirability  function  approach,  desirability  values  close  to  1 are  desired.  The  overall 
desirability,  D",  given  in  Table  3. 10  is  the  geometric  mean  of  the  individual  desirability 
values,  d',  richned  in  equations  (2.4.2)  and  (2.4.3).  In  this  example,  Design-Expert 
yielded  a satisfactory  value  of  0.7018  for  D'.  The  results  for  the  GRG  algorithm 
approach  are  computed  using  the  'solver'  tool  available  in  the  software.  Excel.  For 
details  on  how  to  use  the  solver  tool  of  Microsoft  Excel,  see  Dodge  el  al.  (1993). 
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TVbIt  3.10.  Results  For  the  Desirability  Function  Approach 
RespoDK  SimuH&oeous  Opllma 

1.92  0.52 
19.13  0.94 

yj,  17.57  0.78 

V.,  48.98  0.64 

lootion  » = (-0.6180, 1-000,0.9999)' 

D-  0 7018 


Tabl«  3.11.  FUaulta  for  the  CRG  Approadi 


Response 

Primary  Response 

Location 

1.9177 

1-0.5243,1.4080,0.8618)’ 

20.7467 

(-0.2335.0.4916,0.5218)' 

18.4500 

(-0.4578,1.3580, 0.9728)' 

V4. 

51.2443 

(-0.5842,1.3308, 0.9422)' 

The  results  given  m Table  3.11  are  the  optima  for  each  of  the  four  rcsponeea.  that 
is,  each  of  the  four  responses  was  considered  to  be  the  primary  response  and  the 
optimum  value  was  obulned  using  the  constraints  that  the  other  three  responses  are 
equal  to  their  respective  individual  optima  For  example,  if  jn  is  considered  to  be  the 
primarj’  response,  then  the  maximum  of  yi  is  equal  to  1.9177.  This  optimum  value  is 
obtained  subject  to  the  constraints  that  yu,  and  Vac  given  in  equation  (3.4.13)  arc 
equal  to  18.8366. 17.5833,  53.4116,  respectively-  The  constraint  < 3 was  also 
used.  The  same  procedure  was  foilovred  for  the  other  three  cases,  with  the  constraints 
changed,  depending  on  whidi  variable  is  considered  as  the  primary  response. 

The  optima  using  the  generalised  distance  approach,  the  desirability  function  ap- 
proach and  the  GRG  algorithm  approach,  given  in  Tables  3.9,  3.10  and  3.11  respec- 
tively, are  quite  close  to  each  other.  For  pi,  the  maximum  value  using  the  generalized 
distance  approach,  1.9166,  is  very  dose  to  that  obtained  using  the  desirability  func- 
tion approach,  that  is.  1.9177.  These  two  values  are  also  very  close  to  1.92,  the  value 
obtained  using  the  GRG  appruach.  with  su  , ensideted  as  the  primary  response.  The 


Table  ^^^Th^Ohgina^n^^ode^^evel^Qh^ngi^^^ables 

Coded  Levels 

Variable ^2 ;! 0 I ^ 

A',  65.000  7D.000  75.000  80  000  85.000 

<Vj  4.000  5.000  6.000  7.000  8.000 

.Yj  -0.025  0.075  0.175  0.275  0.375 

A",  0.000  0.0125  0.025  .0375  0.050 

-Vs  0.000  0.050  0.100  0.150  0.200 

Source;  Khuri  and  Cornell  (1996.  p.278) 

tniDimum  values  for  yj  obtained  using  these  three  approaches  are  also  not  that  differ- 
ent. The  same  can  be  said  about  the  minimum  values  of  pg  and  the  maximum  values  of 
y,.  One  should  note,  however,  the  differences  Id  the  way  these  optima  were  obtained. 
-4s  mentioned  In  Ch^ter  2,  Section  2.4,  the  desirability  function  approach  and  the 
GRG  algorithm  approach  do  not  take  Into  consideration  the  correlations  among  the 
responses  and  the  magnitudes  of  tiieir  variances.  These  are  importaoL  considerations 
in  multiresponse  opticmzatlon  (sec  Khuri  and  Cornell  (1996),  p.2SS).  The  generalized 
distance  approach  takes  these  into  account.  Hence  it  is  recommended  that  the  results 
for  the  generalized  distance  approach  be  used. 

Example  2:  Richert  et  al.  (1974)  investigated  the  effects  of  heating  tempera- 
ture (zi).  pH  level  (xs),  redox  potential  (zg).  sodium  oxalate  (zc),  aud  sodium  lauryl 
sulfate  (zg)  on  foaming  properties  of  whey  protein  concentrates.  The  responses  were 
whippiog  time  (pi),  maximum  overruo  (pg),  and  soluble  protein  (pi).  A central  com- 
posite design  consisting  of  a 2^*'  fracliooal  factorial  design,  ten  axial  points,  and  6ve 
center  point  replications  was  lued.  The  results  of  this  study  vrere  utilized  to  And  the 
simultaoeous  maxima  of  all  three  responses.  The  original  and  coded  levels  of  the  five 
ioput  variables  are  given  Id  Table 3.12  and  the  dataset  is  given  in  Table  3.13. 

To  determine  the  best'  multiresponse  modri  in  terms  of  the  PRESS  criterion  and 
fl*.  we  use  the  same  procedure  as  described  in  Example  1.  The  subset  models  given 
Id  Tables  3-14.  3.15.  and  3.16  were  chosen  to  be  the  'best'  (based  on  the  Cp-criterion, 
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Table  3.13.  Experimiintal  Desigp  and  Response  Values  (Example  2) 


Ji  Jj  Xi  Z\  Vi 

in 

y. 

0 0 0 0 0 3.0 

1179 

104.0 

0 0 0 0 0 3.5 

1163 

107  0 

0 0 0 0 0 4.0 

1120 

104.0 

0 0 0 0 0 3.S 

IISO 

101.0 

0 0 0 0 0 3.0 

1195 

103.0 

-1  -1  -1  -I  1 4.75 

1082 

81.4 

1 -1  -1  -1  -1  4.0 

824 

69.6 

-1  1 -1  -1  -1  5.0 

953 

105.0 

11-1-11  0-5 

759 

81.3 

-1  -1  1 -1  -1  4.0 

1163 

80.8 

1 -1  1 -1  I 5.0 

839 

76-3 

-111-11  3.0 

1343 

103.0 

1 1 1-1-1  7.0 

736 

76.9 

-1  -1  -1  1 -1  5.26 

1027 

67.2 

1-1-1  1 1 5.0 

836 

74.0 

-11-111  3.0 

1272 

96.5 

1 1-1  1-1  6.5 

825 

94.1 

-1-1111  3.25 

1363 

95.9 

1-1  1 1-1  5.0 

855 

76.8 

-1  1 1 1 -1  2-75 

1284 

100.0 

1 1 1 1 1 5.0 

851 

104.0 

-2  0 0 0 0 3.75 

1283 

100.0 

2 0 0 0 0 11.0 

651 

50.5 

0 -2  0 0 0 4.5 

1217 

71.2 

0 2 0 0 0 4.0 

982 

101-0 

0 0 -2  0 0 5.0 

664 

85.8 

0 0 2 0 0 3.75 

1147 

103.0 

0 0 0 -2  0 3.75 

1081 

104.0 

0 0 0 2 0 4.75 

1036 

89.4 

0 0 0 0 -2  4.0 

1213 

105.0 

0 0 0 0 2 3.5 

1103 

113.0 

bource;  Khun  and  Uoniell  (1905. 

p.379) 

HifC”’ 

l3l5,I<IS 


3.4.1S  0-93827  0-31482  21-000  i|,ra,ii,i4,  xs, 

3,.,..  ppp- 


sas:"'" 


Tsbie  3-l^Mul^lgg^^ModeU_^mc^^sin^Mod^3AH  to  3.4.25 


No.  for  y,  for  U2  foe  yj 

1 3.4.14  3.4.1S  3.4.22 

2 3.4.15  3.4.18  3.4.22 

3 3.4.16  3.4.18  3.4.22 

4 3.4.17  3.4.18  3.4.22 
3 3.4.14  3.4.19  3.4.22 

6 3.4.15  3.4.19  3.4.22 

7 3.4.16  3.4.19  3.4.22 

8 3.4.17  3.4.19  3.4.22 

9 3.4.14  3.4.20  3.4.22 

10  3.4.15  3.4.20  3.4.22 

11  3.4.16  3.4.20  3.4.22 

12  3.4.17  3.4.20  3.4.22 

13  3.4.14  3.4.21  3.4.22 

14  3.4.1S  3.4.21  3.4.22 

15  3.4.16  3.4.21  3.4.22 

16  3.4.17  3.4.21  3.4.22 

17  3.4.14  3.4.18  3.4.23 

18  3.4.15  3.4.18  3.4.23 

19  3.4.16  3.4.18  3.4.23 

20  3.4.17  3.4.18  3.4.23 

21  3.4.14  3.4.19  3.4.23 

22  3.4.15  3.4.19  3.4.23 

23  3.4.16  3.4.19  3.4.23 

24  3.4.17  3.4.19  3.4.23 

25  3.4.14  3.4.20  3.4.23 

26  3.4.15  3.4.20  3.4.23 

27  3.4.16  3.4.20  3.4.23 

28  3.4.17  3.4.20  3.4.23 

29  3.4.14  3.4.21  3.4.23 

30  3.4.15  3.4.21  3.4.23 

31  3.4.16  3.4.21  3.4.23 

32  3.4.17  3.4.21  3.4.23 


2154.3046 

2027.3294 


1947.4140 

1803.3604 

1301.7832 

1433.4466 

1616.1234 

1513.9965 

1478.3295 

1037.0494 

1610.0459 

1506.0677 

1471.9379 

1021.4802 

1885.5151 

1751.3511 

1749.9175 

1390.0528 

1604.4202 

1455.3973 

1454.0606 

1068.8917 

1259.5952 

1149.4299 

1116.0802 

640.1738 

1257.3202 

1145.3039 

1113.5910 

630.3612 


0.9445 

0.9433 

0.9433 

0.9246 

0.9250 

0.9421 

0.9421 

0.9220 

0.9225 

0.9223 

0.9228 

0.8996 

0.9225 

0-9197 

0.9202 

0.8970 

0.9430 

0.9417 

0.9417 

0.9206 

0.9223 

0.9405 

0.9405 

0.9178 


0-9167 

0.9172 

0.8903 
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Mod«l 

■33 

34 

35 

36 

37 
3S 
39 

41 

42 


45 

46 

47 

48 

49 

50 

51 

52 

53 

55 

56 

57 


60 

61 

64 


Model  Modri  Model  PF 

tor  ai  for  !n  to  W 


3.4.14  3.4.18 

3.4.15  3.4.18 

3.4.16  3.4.18 

3.4.17  3.4.18 

3.4.14  3.4.19 

3.4.15  3.4.19 

3.4.16  3.4.19 

3.4.17  3.4.19 

3.4.14  3.4.20 

3.4.15  3.4.20 

3.4.16  3.4.20 

3.4.17  3.4.20 

3.4.14  3.4.21 

3.4.15  3.4.21 

3.4.16  3.4.21 

3.4.17  3.4.21 

3.4.14  3.4.18 

3.4.15  3.4.18 

3.4.16  3.4.18 

3.4.17  3.4.18 

3.4.14  3.4.19 

3.4.15  3.4.19 

3.4.16  3.4.19 

3.4.17  3.4.19 

3.4.14  3.4.20 

3.4.15  3.4.20 

3.4.16  3.4.20 

3.4.17  3.4.20 

3.4.14  3.4.21 

3.4.15  3.4.21 

3.4.16  3.4.21 

3.4.17  3.4.21 


3.4.24  2050.8224 
3.4.24  1913.5294 
3.4.24  1912.5955 
3.4.24  1583.4237 
3.4.24  1812.6690 
3.4.24  1660.0629 
3.4.24  1059.2509 
3.4.24  1307.8545 
3.4.24  1489.7999 
3.4.24  1382.8109 
3.4.24  1348.1879 
3.4.24  911.1950 
3.4.24  1485.1354 
3.4.24  1375.5986 
3.4.24  1342.6031 

3.4.24  S97.29S4 

3.4.25  1622.0860 
3.4.25  1506.4869 
3.4.25  1507.8704 
3.4.25  1037.6894 
3.4.25  1269.4878 
3.4.25  1143.1515 
3.4.25  1144.8639 
3.4.25  688.5999 
3.4.25  950.4228 
3.4.25  847.2803 
3.4.25  827.7270 
3.4.25  256.2547 
3.4.25  951.0136 
3.4.25  845.1851 
3-4.25  826.8827 
3.4.25  249.5952 


0 9433 
0.9417 
0.9417 
0.9213 
0.9225 
0.9405 
0.9405 


0.9198 

0.9193 

0.9198 


0,9171 

0.3906 

0.9275 

0.9258 

0.925S 


0.9242 


0.8554 


0.8517 


Table  3ii8i  IhgSaMtis 

0.S3375  -12.65426  -2.45327 

-12.65426  5048.52800  23.15932 
-2.45327  23.15932  38.71550 


ble  3.18.  FVcm  these  values  is  computed  yielding 


= (3.9474, 1.2708, 0.1875,  -0.4375,  -0.1875,  -0.0833, 

0.8805, 0.7191,  -0.3361,  -0.3277)', 

= (1115.396,-176.0833.-18.1667,57.5833,21.8333,19.0833, 
-50.4371,  -37.5328,  -44.5806)', 

■y,  = (101.8229,  -8.2458, 7.5125, 2.3792. 1.1292. 1.6625, 

-7.2013,-4.1430,-1.5686)'. 


Using  the  7^  values  we  have 

iu  = 3.9474-H.2708H -•■0.1875*3-0.437513  - 0.18751, -0.0833X4 -I- 
0.88051?  -I-  0.7191x,i3  - 0.3361*3ij  - 0.3277i,is 
ya,  = 1115.396-  176.0833*1 -18.16671, -•■57.5833x3  + 21.83331, -H9.0833xs- 
50.4371X?  - 37.5328*?  - 44.5806*1*3 

S3.  = 101.8229 -8.2458*, + 7.5125*3  + 2.3792*3 + 1.1292*. + 1.6625*3  - 

7.2013*?  - 4.1430^  - 1.5686*?.  (3.4.28) 

The  individual  maxima  are  again  obtained  using  the  FORTRAN  program  by 
Conion  (1992)  (see  Appendix  B,  Section  6.3).  These  maximum  values  and  their 
locations  are  ^vea  In  Table  3.19.  The  experimental  region  is  given  by  E?,,  i?  < 5. 
The  minimum  of  the  disuoce  measure  p given  in  equatinn  (3.3.1)  is  found  based  on 
26,000  values  nf  p using  the  procedure  PROC  IML  of  SAS.  Its  value  aod  locatioo  are 
giveo  in  Ihble  3.20.  The  simuluneous  maxima  are  also  given  in  Table  3.20.  These 

•vith  the 


than  their  I 
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Table  3.19.  The  Individual  Maxima  of  the  Three  Responses 

Raponae  Individual  Maxima Location 

])„  «,  = 11.8031  (1.9991.0.9248.-0.3550.-0.1378.0.0093)’ 

yj.  4,  = 1385.7395  (1,4292.0.3501,1.1633.0.9862,0.1533)’ 

e,  = 113.5343  [-0.5185,0.7689,1.6597.0.2509.1.1498)' 


_^ble^jffij_TheJimulUBeoua^l^in^^ 

Reipoiiee Simultaneoua  Maxima 

5.9897 

,j„  1289.1471 

78-9504 

Location  a = (-1.9103, -1.0731, -0.0287, -0.0390.0.3628V 
minp 17.2141 


to  78.9504.  The  minimax  approach  described  in  the  previous  section  can  also  be  used 


CHAPTER 


MULTIRESPONSE  MODELS  WITH  FLXED  BLOCK  EFFECTS 
4.1  Introduction 

Far  the  most  part,  analyses  results  available  In  the  statistics  litetature  conceming 
response  surface  models  with  blocli  effects  apply  only  to  the  single-cespouse  case.  In 
experimental  situations  where  several  responses  of  interest  are  measured  and  where  it 
is  impossible  to  perform  all  the  runs  under  homogeneous  conditions,  the  analysis  of  the 
resulting  multiresponse  data  must  take  into  account  heterogeneity  of  the  experimental 
conditions.  For  example,  a diemic^  engineer  may  choose  to  run  an  experiment  using 
two  separate  machines  (instead  of  only  one)  located  in  tura  different  buildings.  In 
this  case,  machine  confounded  with  location  should  be  included  in  the  analysis  of  the 
experiment.  This  extraneous  source  of  variatiOD,  that  is,  the  effect  of  location  (block 
elfecti,  can  Influence  the  results  of  the  analysis  of  the  data,  specifically,  the  estimation 
of  the  mean  response  over  a given  region  of  interest.  Thus  such  block  effects  must  be 
considered  in  any  multiresponse  analysis  of  the  data. 

In  this  chapter,  the  results  ^ven  by  Khuri  (L994)  coDceming  the  effect  of  blocking 
on  the  estimation  of  a response  surface  are  extended  to  the  multiresponse  esse.  Here 
we  consider  the  estimation  of  the  model  parameters  aod  of  the  mean  response,  the 
effect  of  blocking  on  the  pwidiction  variance,  and  present  a condition  for  orthogonal 
blocking.  These  are  all  important  considerations  in  the  analysis  of  rssponsc  surface 


experiments  with  block  effects.  The  model  assumptions 
(1994)  and  these  are  given  in  the  next  section. 


follow  Khuri 


Khuri  (1994)  discussed  the  effect  of  blocking  on  the  estimation  of  a response 
surface  for  the  single'tesponse  case.  He  considered  the  model 

y-  = JaU  + X-0-  -t-Z’r  +e-  (4.2.1) 

wliere  y'  is  the  response  vector  which  consists  of  n observations;  1,  la  a vector  of  ones 
of  order  ti  x 1:  X'  is  an  n k p matrix  whose  uth  row  is  where  x„  is  the  vector 

of  design  settings  at  the  uth  experimental  run  (u  = 1, 2,. . . » a)  and  / (Xi.)  is  a vector 

of  order  Ixp  whose  elements  consist  of  thex.i  terms  ((  = 1,2 k)  along  with  their 

powers  and  cross*products  of  these  powers;  is  a vector  of  unknown  itarsmetets; 

6'  ^ {6\,6i de)',  where  Sj  denotes  the  effect  of  the  jth  block  ( j=l,2,  ...  ,b);  and 

Z'  is  i block'diagonal  matrix  of  the  form  Z = dio^(lAt,  Iq Ing),  where  Uj  Is 

the  Else  of  the  jth  block,  that  is,  the  blocks  are  of  different  sizes.  The  random  error 
vector  e'  is  assumed  to  have  a zero  mean  and  a variance-covariance  matrix  a‘I„. 

Model  (4.2.1)  can  be  written  in  a more  compact  form.  Since  the  columns  of  Z' 
sum  to  In,  model  (4.2.1)  can  be  written  as 

p-aX’/r+ZV  + e-  (4.2.2) 

where  t’  = Al,  + 6'.  If  the  block  effects  are  constrained  to  sum  to  zero,  then 
tk  = {l/b)  E5=i  , where  t'  is  the  jth  element  of  t'.  .Model  (4.2.2)  can  be  written  as 
y-  = W 9‘  + (•  (4.2.3) 

where  W = (X' : Z’lnnd  8*  )'. 

We  now  extend  the  above  results  to  the  multiresponse  case,  that  is  obtain  the 


el  and  determine  the  estimates  of  the  model's  parameters.  Suppose 
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there  are  r responses  of  interest  and  that  all  the  reponscs  have  the  same  number  of 
observations,  ti.  The  polynomial  models  for  the  different  responses  can  have  different 
degrees,  d,.  In  some  tnuitiresponse  experiments  it  is  possible  that  the  blocks  used  are 
not  necessarily  the  same  for  all  the  responses.  For  example,  a chemical  engineer  in  a 
production  plant  wants  to  improve  the  operating  conditions  for  the  production  of  a 
certain  polymer,  that  is,  he/she  warns  to  determine  conditions  on  the  input  variables 
temperature,  pressure,  conceutrallon  of  material  A,  and  dye  concentration  that  will 
optimize  the  responses  given  by  mechanical  strength  and  coloration  of  the  polymer. 
.Note  that  material  A is  used  to  affect  the  mechanical  strength  of  the  polymer  and 
does  not  really  affect  coloration  and  that  the  dye  does  not  really  affect  mechanical 
strength.  The  engineer  found  out  that  there  are  two  batches  of  material  A and  three 
batches  of  dye  available  for  the  experiment.  In  this  case,  two  different  blocks  may 
be  used  for  the  two  responses,  that  is,  the  2 batches  of  material  A for  mechanical 
strength  and  the  3 batches  of  the  dye  for  coloration.  For  example,  if  a CCD  for  it  = 4 
factors  is  used  to  run  this  experiment,  then  for  mechanical  strength,  the  two  blocks 
given  by  the  factorial  portion  with  center  point  replicates  and  the  axial  portion  with 
center  point  replirates  may  be  used.  For  coloration,  two  of  the  three  biocks  each 
consists  of  a half  replicate  of  the  2^  factorial  points  plus  center  point  replicates  and 
the  third  block  consists  of  the  axial  points  plus  center  point  replicates.  Thus  in  setting 
up  the  multirespoose  model,  we  will  consider  the  case  in  which  the  blocks  are  not 
necessarily  the  same  fur  all  the  responses.  From  models  (4.2.1)  and  (4.2.3),  the  Ith 
single-response  model  is  given  by 

y,  = iiotln  * X,0,  + Z^6i  + 1 = 1 r,  (4.2.4) 


equivalently, 
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whew  fVi  = [Xj : Z,|;  Z,  = <(toj(U„,l„„ In,,,),  where  n,j  is  ihesize  of  the;th 

block  for  the  ith  response.  5,  is  tire  number  of  blocks  for  the  ith  response,  Sjsi  ~ 't' 
I = 1.2 r;  X,  is  of  order  n » p,;  9,  = {(?,',  0,  is  » vector  of  unknown 


of  the  jlh  block  in  the  ilh  response  model,  j = 1,2,  — 6,  -,i  = 1,2 r.  It  is 

assumed  that 

£(e,)  = 0,  i = l,2 r,  ond  £(e,ei)  = <r,i/„,  ».i  = l,2 r.  (4.2.6) 

The  response  models  in  (4.2.5)  can  be  expressed  as  a single  multireaponse  model 
of  the  form 

y=H'8  + e (4.2.7) 

vrhere  y = (pi,  V}, . (V  is  a block-diagonal  matrix  of  the  form 

IV  = diagiWx.Wi W,).  9 = and  e = (e',,ei <)'-  By 

(4.2.6),  the  variance-covariance  matrix  of  e is  where  £ ^ (nry)  is  the  variance- 

covariance  matrb!  of  the  r responses. 

The  estimates  of  the  parameters  in  model  (4.2.7)  can  be  obtained  by  using  results 
from  Srivastava  and  Giles  (1987,  p.12).  The  generalized  least  squares  estimator  of  9 


8 = (IV'(£-'  ®/„)H')-'H'‘(£-'  ®/„)y  (4.2.8) 

and  its  vanance-covariance  matrix  is 

b'or(81  = l»"(E-'®/„)»’|-'.  (4.2.9) 

Let  E"‘a/„  = (o''/,).  Then,  using  the  result  given  by  Zellner  (1962,  p.351,  formula 
2-8),  the  variance-covariance  matrix  of  8 can  be  written  in  the  form 


l'ar(8) = 


<r“iv;»'i  rr'^H'^H',  ...  a''W\W,  y' 

.■•"W'iW,  n”W'.,Wi  ...  a‘’W'jW, 

a'>W,W,  n'^W'.W,  ...  j 


(4.2.10) 


1q  geoeral.  the  variance-covariance  matrix  T,  is  unknown.  To  obtain  the  estimates 
in  (4.2.8),  S must  be  estimated.  To  estimate  £.  we  ran  use  the  estimator  given 
Zellner  (1962)  (see  also  Srivastava  and  Giles  (19S7,  pp.l2-i7)).  This  eslUnator  was 
denoted  in  Chapter  3 by  5 = (Sg)  where 

»u  = - H',(tv;iv,)-'H';ii/-  iv,(w;tVj)'‘»’;iv,-.  (4.2.11) 

Using  5 in  place  of  £ in  (4.2.8),  an  estimated  generalised  least  squares  (ECUS) 
estimator  of  $ is  given  by 

8,  s (H"(S-'  ®/„)1Vj-'W'’(S-'®/„)y.  (4.2.12) 

This  estimator  does  not  have  the  same  properties  as  8.  Unlike  8 in  (4.2.8),  6,  no 
longer  has  the  minimum  variance  property.  From  (4.2.9),  an  approximate  expression 
For  the  estimated  variance-c  "variance  matrix  of  8.  is  given  by 

l'ir(8,)  = (Wr'(5-i  »/„)»■)-'.  (4.2.13) 

Flirthermore,  if  X,  - X,  and  Z,  = which  means  that  Wi  = W,  for  i = 

1.2 r,  then  model  (4.2.7)  takes  the  form 

y = {l,eiW,)e  + e.  (4.2.14) 

In  this  case,  the  CLS  estimator  of  8 is  ^ven  by 

8 = {/r  ® |(K'’.»'.)-'lV;]}y  (4.2.15) 

which  is  the  form  of  the  ordinary  least  squares  (OLS)  estimator  of  8. 

4.3  Estimating  the  Mean  Response  in  the  Pr>«>tv.p  of  Fixed  Block  Effects 

In  Section  3.1,  estimators  of  the  model  parameters  was  obtained  by  using  the  GLS 
method.  In  this  section,  we  derive  the  estimators  oi  the  mean  responses  by  applying 
the  same  approach  and  by  using  the  results  in  (4.2.8),  (4.2.9),  and  (4.2.12). 
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Suppose  lhal  the  ilh  response  surface  model  is  a polynomial  of  order  d,  of  the 

<il*)  = (4-3.1) 

where  2 s {x,,X2, . . . ,x«)'  is  a vector  of  k input  variables.  t7i(x)  denotes  the  mean 
of  the  ith  response  at  x,  and  f',{x)  is  a vector  vrhose  elements  consists  of  Che  z,'s 
defining  the  polynomial  temui  up  to  degree  dt.  The  elements  of  and  dot  sre  unknown 
parameters.  The  least-squares  estimator.  ^,(x).  of  the  Ith  response  at  x is 

^(x)  = jci,  + /|(x)ft  (4.3.2) 

where  ^ Based  on  this  result,  formula  (4.3.2)  can  be  written  aa 

ht(x|  = 9i’(x)e,  (4.3.3) 

where  is  the  unbiased  estimator  of  8,  obtained  from  (4.2.8)  and  pi'(x)  = (/r'(x) : 

plfcl-  ^(®)  = (i7i(x),%(x) flrix})’.  Using  (4.2.9),  the  variance-covariance 

matrix  of  i^{x)  is 

Vnr(^(*)]  = ^[^"(E-i  a I„)W]-‘G  (4-3.4) 

where  G = dioj(9,(x),  j.;(x) g,(x)).  If  £ is  unknown,  we  use  the  EGLS  estima- 

tor. 6„  given  in  equation  (4.2.12),  that  is, 

%{x)  = C e„  (4.3.5) 


‘'“'•|n.(x)J  = G-[»”(5-'  a I„)W]-'G. 


(4.3.6) 


Khuri  (1994)  pointed  out  in  the  single-response  case  that  the  block  effects  cause 
an  increase  in  the  prediction  variance,  assumiog  that  the  variance  of  the  experimental 
error  term  in  the  model  is  the  same  whether  or  Dot  block  effects  ace  Included  In  the 
model,  under  the  condition  that  blade  effects  potentially  exist.  Specifically,  if  block 
effects  exist,  it  Is  assumed  here  that  the  variance  of  the  experimental  error  term  is 
the  model  with  blocking  terms  is  the  same  as  the  variance  of  the  experimental  error 
term  in  the  model  without  blocking  terms.  Khuri  (1994)  shovred  that  for  a particular 
response,  say  response  i, 

Vor(^(x)|  > /ar[pi,(x)| 

where  4(z)  follows  the  form  given  in  (4.3.3)  and  ^(x)  is  the  least-squares  estimator 
of  the  ith  response  at  x when  the  block  effects  are  zero.  He  proved  this  inequality  by 
establishing  the  identity  given  by 

I'orl^li))  = Vor(^(x)|  -+-  ff;(x)Q‘g,(x)er’.  (4.4.1) 

In  the  above  equation,  Var[^(x)|  follows  the  form  co’,  where  c = ft;(x](t/|y,)‘‘fi,(x), 
hKa)  = [I  : V.  = |1„  : X,],  g((x|  is  given  in  (4.3.3),  nU,  is  the  variance- 

covariance  matrix  of  si  io  (4.2.4).  Q’  is  giveo  by 

Q-  = lW,W,]-'K.\K',{W',W,y'K.]-'K‘,[W.W,)-'.  (4.4.2) 

W,  is  defined  In  (4.2.5),  w [0  : L,\,  it  = |1,._, ; and  0 is  a zero  matrix 

of  order  (5,  - 1)  x p,,  with  p,  being  the  number  of  columns  of  X,.  This  Implies  that  a 
model  containing  blocking  terms  will  have  a prediction  variance  value  that  Is  greater 
than  or  equal  to  the  prediction  variance  for  a model  without  blocking  terms  since  Q' 

Is  positive  semidefipite.  Thus  the  addition  of  the  blocking  terms  will  not  reduce  the 
prediction  variance,  assuming  that  the  variance  of  the  error  term  for  the  blocked  data 
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is  the  ssme  as  the  variance  of  the  error  term  for  the  unbiocked  data.  Id  this  sectioo. 
we  deuioDstrate  that  the  same  is  true  for  the  muitiresponse  case. 

Let  j = (j',.  jj,  ....f,)',  where  5,  is  the  vector  of  block  effects  for  the  ith  respoose. 
Let  Vorli^gfx))  be  the  predictioo  variance  when  the  block  effects  are  zero,  that  is, 
the  eiemeots  of  S are  zero.  We  will  show  that 

V'or(n(x)|  > Var[T|o(®)) 

for  all  s in  a re^on  of  interest  !R.  that  is,  I^ar[^(x)|- V’ar|q,(z)]  iapoeitivesemideli- 

niteforallz  € ?!.  To  show  this,  we  note  that  if  j = 0,  then  r,  = i = 1,2 r. 

This  means  that  the  elements  of  Ti  (i  = 1,2,, . . , r)  are  equal.  We  can  therefore  write 
tiT,  = 0 where  L,  = [l^-i  : -/(,-iJ.  Also,  the  prediction  equation  for  the  ith  re- 
sponse can  be  expressed  as  ^(x)  = gi{x)  where  9^^  is  the  least  squares  estimator 
of  fit  under  the  restriction  L,r,  = 0 or  equivalently  K'0,  s 0 where  K\  is  defined 

in  (4.4.2).  Hence,  K'8  = 0 where  K‘  =<HaglK\.K'2 K',).  Now,  iisin|  a result 

given  by  5eber[1954,  formula  8.23,  p.403),  the  generalized  least  squares  estimator  of 
6 under  the  above  restriction  is 

$0  = 9 - A-'K(K‘A  'K)-'K‘e  (4.4.3) 

where  A = IV'fE"  ® In)W.  From  this  result,  the  variance-covariance  matrix  of 
00  is 

Var(9o)  = PA'P-  (4.4.4) 

wliere 

P = I-A-'KiK‘A'K)-'K' 

This  variance-covariance  matrix  can  also  be  written  in  the  form 

i^ar(9o)  = A-'-A-'KtK‘A-'K)-'K'A  '-A-‘K(K'A''K)  K’A  ' 
■t-A-'K{K  A 'Ky'K'A-'KlK'A-'Ky'K'A-^ 


= A-'  - A-'K{K‘A-'K]-'K‘A-'  - A-'KIK'A-'K)-'K'A-' 
+A-'K{K'A  'K]  'K'A  ' 

= A-'  - A-'K{K‘A-'K)-'K'A-'. 

The  variance-rovariDDce  matrix  of  17(3;)  in  (4.3.4)  can  then  be  expressed  as 
i-or[«x)j  = GA-'G 


where  Q = A~' K{K' A~' K)~' K' A~'  is  of  order  (p  + n=i  x (P+  P = 

pi,  and  the  rank  of  Q is  equal  to  the  rank  of  K.  which  is  equai  to  bi-r. 

Since  Q is  positive  semidefinitc.  it  foliows  that  Varj^fx))  > Var|q,(z)|,  assuming 
that  S remains  unchanged.  However,  if  E is  replaced  by  its  estimator  S,  then  this 
inequality  may  not  necessarily  hold. 


It  is  often  the  case  that  the  polynomial  effecis  In  the  models  described  in  equa- 
tion (4.2.4)  are  not  independent  of  the  block  effects.  Box  and  Hunter  (1957)  pointed 


blocks  greatly  affects  the  efficiency  of  the  estimation  of  the  model  parameters.  .4s  a 
result,  it  would  be  desirable  to  analyse  the  polynomial  effects  independeotlv  of  the 
block  effecis.  We  can  achieve  this  by  ortJiogmat  Hoektnp.  A design  in  which  the 
polynomial  effects  can  be  estimated  independently  of  the  block  effects  is  said  to  block 
orthogonally.  Thus  Knding  conditions  for  orthogonal  blodting  for  the  multireaponse 
case  is  desired.  Consider  the  models  given  in  equation  (4.2.4).  The  multiresponse 
form  of  these  modeb  is 


= G lVor(9o)  + 0]  G 
= Vor[qo(x))  + G'OG 


out  that  the  manner  in  which  tlie  experimental  units  of  a design  are  assigned  to  the 


V = Mfti  + X3  + 2i 


(4.5.1) 


M = 1,0  U.  00 


where  Aa  = dft  + and  Z.  = (/„  - 1 J„}2.  The  muliiresponse  form  of  the 


y = MAo  + XB  + Zd  + 


= (M  ; A-)  and  0‘  = (A,' : 0)'. 


of  j.  On  thIebE 


X'^(Z-'  0l„)2  = 0. 

et  E-'  = (ff-J),  i.j  = 1,2 r.  Then  equation  (4.5.5 

^ S:l: :::  S;!:! 


c.‘'x;z,  n'^X'.Z-^  ...  a'’X\Z. 

c"X':Z, 

0''x;z,  o’^X'.Z:  ...  n"X,Z.) 
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Thus  If  condition  (4.S.5)  or  14.5.6)  issstisficd,  it  is  possible  to  estimate  the  polynomial 
effects  in  inultirespntise  model  (4.5.4)  independently  of  the  block  effects.  Based  on 
this  condition  for  orthogonal  blocking,  two  theorems,  one  on  the  estimation  of  model 
parameters  and  the  other  on  the  prediction  variance  will  be  established.  The  frdlowing 
results  are  needed  to  prove  these  theorems;  from  equation  (4.5.6),  the  matrix 

/ tr“i;,(J,/n)Z,  o'»i;.(J„/n]Z, 

<^'V„[J^ln)2,  o“i;,(J„/n]Z, 

o''i;,(J„/n)Z,  o'»l',(J./n)Z, 

o"X;(J,/n|Z,  a'^X\{J,ln)Z, 

<r''Xi(J,/n)Z,  ir“X’j(J,/n)Zj 

, o''X;(J,/n)Z,  a'^X',[JJn)Zt 

is  equal  to 

®/„ 

X’(E-'0/„ 

Let  J = I,  % J„.  From  (4.5.7)  and  (4.5.8)  we 
X'(E-'®/„)Z  = X'(E 

If  condition  (4.5.5)  is  satisfled.  some  important  results  can  be  obtained.  To  be 
able  to  derive  these  results,  we  hist  consider  Lire  models  described  in  equation  (4.2.5), 

B,=  X,A  + Z,t.4-€„  i=l,2 r 

where  t,  = Hi,!,,  + S,.  The  multiresponse  form  of  the  above  modris  is  then 

p = X(J  + Zt  4-e  (4.5.10) 

where  r = (t'j.Tj r^)'  and  X./3.Z,  and  e are  as  defined  In  model  (4.5.1).  Note 

that  the  latter  model  is  the  same  as 


o^'i;(J,/n)Z, 

o‘'X'(Jn/n)^, 

rr^X',(J„/n)Z, 

ff"X;(J„/n)Z,  ) 


= X'd  + Z«  + e 


(4-5.11) 


where  X*  is  as  defiDed  io  equatioo  (4.5.4)  end  i9  s : ff)'. 

As  stated  earlier,  urlliogonal  bluckiog  allows  us  to  estimate  the  polyoomial  effects 
iodependemly  of  the  block  effects.  In  the  following  theorem,  we  also  show  that 
orthogonal  blocking  allows  us  to  get  correct  estimates  for  the  polynomial  effects  even 
if  the  block  effects  were  incorrectly  excluded  in  the  analysis  of  the  multiresponse 
experiment. 

Theorem  4.J.J  //  (he  design  6focjks  orthogoruiUyr  (hen  the  genemheed  least  aquares 
eattmalor  0 of  0 in  modcf  (4-1.10)  is  (he  same  as  when  the  block  effects  S are  lero. 
that  IS.  0 = 0 isArrr  0 satssjies 

^ j = |X''(E-‘®I„)Xn-'X^(E-'  ®/.)W.  (4.5.12) 

Proof;  The  proof  is  given  in  Appendix  0,  Section  D.l. 

The  importance  of  the  above  theorem  lies  in  the  fact  that  even  if  an  experimenter 
incorrectly  analyses  a multiresponse  experiment  with  block  effects  by  excluding  them 
in  the  analysis,  the  resulting  CL.S  estimator  for  0 would  still  be  correct,  provided 
that  the  design  used  blocks  orthogonally.  The  other  important  result  for  orthogonal 
blocking  is  established  In  Theorem  4.5.2.  In  Section  4 it  vras  pointed  out  that  the 
block  effects  cause  an  increase  in  the  prediction  variance.  In  the  following  theorem 
we  show  that  when  the  design  blocks  orthogonally  the  prediction  variance  at  a point 
X Inside  the  experimental  region  exceeds  l^ar[^,(zj]  by  a constant  amount  which 
depends  on  the  sizes  of  the  blocks  for  a given  E. 

TTieorrm  2.5.3  If  the  design  blocks  orthogonally  and  if  the  blocks  are  the  same /or  off 
the  T responses  of  tnlerasl.  (fial  is.  = . . . = Z,  = Ze,  then  the  pirdiclion 
variance  jmen  m eguation  I4.S.4)  takes  the  form 
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•oherr  V'nr[iio(*)]  iJtrtiolpj  the  predtiiioti  variance  luAen  block  effects  d = 0,  Uiat 


where  H‘  = [/. ; fl  f”  = <hag{f\lx),fiix) /U®))'  /!(*)  “ ““  hefinai  \ii 

msiie/  (4-H-l/- 

Proof;  Thf  proof  is  givon  in  Apppndij  D,  Sprtion  D.2. 

Note  that  to  apply  Tlieorem  4.5.2.  we  need  not  aesunie  that  the  polynomial  models 
for  the  responses  hate  the  .same  tlegree.  that  is,  the  polynomial  models  ran  have 
different  degrees  d,. 


In  this  example,  we  use  the  results  of  an  experiment  performed  by  Chen  el  al. 
(1989)  cutirerniiig  the  poet-storage  rendition  of  freshly  harvested  Delicious'  apples. 
The  apples  were  harvested  from  the  same  orrharris  at  tsro  different  dales,  that  is, 
‘Harvest  I'  is  on  day  I of  the  t;ommert;ial  harvesting  date  and  ‘Harvest  IF  is  14 
days  after.  We  consider  these  two  harvesting  dates  as  the  fixed  block  effects  in  the 
experiment.  The  input  v‘ariables  considered  were  Xi  — number  of  weeks  in  storage 
after  harvest  and  A'j  = storage  temperature  (*C)  after  harvest.  There  were  three 
response  variables  of  intertat;  ijt  = Hesh  Hrinness  (N).  pj  = amount  of  extractable 
Juice  (ml/IOOg  fluid  weight),  and  !/]  = internal  ethylene  composition  of  the  apples 
(ppm).  A 4 X 4 factorial  design  was  used.  The  design  settings  in  the  uriginal  and 
coded  variables  and  the  multiresponse  data  are  given  in  Table  4.1.  Note  that  the 
square  root  of  the  m values  are  given  in  Table  4.1  and  that  the  coded  values  of  A', 
and  A'j  are  derived  from  the  following  equations; 


-nrfedi]  = H'|AC'(2-‘  s;/„)X*]-'H, 


(4.5.14) 


■Vi  - 2.5 
1.12 
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OrijcinAl  Input  \'itriAhlf«  Corletl  Input  VarUbln 


-1.21  62.41  75 


-0.64  57.74 
-D.64  58.18 
-0.64  58.89 
-0.64  55.29 
0.45  57.71 
0.45  51.77 
■ 41.77 
0.45  .28.25 

1.38  53.62 

1.38  43.59 

1.38  36.52 

1.38  32.01) 
1.21 

1.21  59.65 

1.21  58.34 

1.21  57,02 


; 68  12.0108 


■Ya  - 9,3 
7 70 


where  2.5  and  9.3  are  the  means  of  -Yi  and  X-,,  respenively,  and  1.12  and  7.7(1  are 
the  rorreeponding  standard  deviations,  where  the  formab  . i = 1,2, 

was  used.  From  Table  4.1  we  ran  see  that  all  the  responses  liave  the  same  blocking 
arrangement. 

For  this  example,  we  wish  in  determine  optimum  operating  conditions  on  .storage 
rinte  and  temperature  for  both  harvest  rimes.  To  do  so,  we  first  need  to  determine 
the  'best'  multiresponse  model,  taking  into  account  the  fixed  biock  effects.  Once  the 
model  Is  established  we  ran  use  the  optimization  procedure  given  in  Chapter  3 to 
obtain  the  simuicajieuus  optima  of  the  three  responses. 

Using  the  Cp,  and  adjusted  criteria  described  in  Chapter  3,  the  ‘best’model 
for  each  response  was  determined.  The  proredore  PROC  RSQUAREofSAS  was  used 
to  detennine  these  iiiodeis  and  the  results  are  given  in  .Appendbt  E.  A square  root 
transformation  was  itsod  to  stabilize  ths  residuals  for  ps.  This  means  tliat  the  vector 
used  in  the  model  and  parameter  estimation  contains  the  square  roots  of  the 
observed  values  for  trj.  .A  t'oinpleie  second-order  model  in  X|  and  xe  was  considered 
to  be  the  best'  for  each  of  the  three  responses  lisaed  on  the  Cp,  and  adjusted  K‘ 
values  given  in  Appendix  D.  .Applying  e(|uation  (4.2.5),  the  'best'  model  for  ail  three 
responses  can  be  written  in  the  form 

V,  ‘ + 1 = 1.2,3 

where  W„  = [X.  : Z„\.  Z,  = dins(li«. lis),  and  X,  Is  a matrix  whose  utb  row  is 

given  by  /‘(x,)  = ii  = 1.2 32.  Using  equation  (4.2.14) 

the  multiresponse  form  of  tliese  models  is  given  by 
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Tatilf  4.2.  Cp.  H‘,  and  adjusted  fl’  values 
Criterion  qi  1)2  1)3 


0.936  0.979  0.930 
adjusted  R*  0.923  0.975  0.916 
Cp 1 G 6 6 


where  y = (vl.  yi,  s^)‘-  9 = (S'!.  and  e = (s',,  e^,  ^3)'.  Note  that  since  Wi  = 

W2  = IV,  = W„  the  GLS  estimator  of  fl  does  not  depend  on  S,  that  is,  the 
GLS  estimator  roinrides  with  tlie  OLS  estimator.  Srivaaiava  and  Giles  (1987,  p.lT, 
etjuations  2.15-2.16)  showed  that  in  this  ease,  the  EGLS  estimator  is  the  same  as  the 
OLS  estimator.  Titus  using  the  OLS  estiitiator,  $,  = for  t = 1,2,3, 


Using  equation  (4.3.3),  tlie  least  squares  estimates  of  qi(x}  are 

ill*)  = 50.09  -4.13i|  - 7.79i,  + 0Jli^-e  l.05i|  - 2.85iiij 
i)j(x)  = n.59  - 0.96i,  - 1.53ij-i-0.24x?-l-0-42^ -O.GliiX, 
ij,(x)  = 8..395  + 0.72i,4-.3.57i,-(l.54i?-0.87*|  + 0.33i,ij. 

Note  tliat  it(x)  here  corresponds  to  the  .square  root  transformation  used  for  ys.  The 
Cp.  !P.  atid  adjusted  R’  for  the  above  models  are  given  in  Table  4-2- 

It  is  desired  to  Hnd  the  operating  conditions  that  tnaximize  y,  and  pj,  atid  minimize 
P3-  For  this  purpose  we  can  apply  the  simultatieous  optimization  procedure  given  in 
Oiapter  1 Tlus  amounts  to  minimizing  the  distance  function  given  by 


(-4.13,-7.70,0-31,1.05,-2.85,50-61.49.57)' 


( -0.96,  -1.53, 0.24, 0.42,  -0.61, 71 .84. 71 .34)’, 


(0.72.3.57,  -0.54,  -0-87, 0.33, 8.30,8.49)’. 


where  i>  is  ihe  r X 1 esiimatefl  iwtor  of  individual  optima  over  the  experimental 
region  given  by  the  inequalities  —1.34  S >t|  $ 1-34,  and  -1.21  < xj  < 1.33, 
g'(®)  is  a vector  of  order  1 y p,  that  has  the  same  form  as  a row  of  W„  but  Is 
evaluated  at  the  point  x.  is  tiie  number  of  columns  of  W^.  and  £ = - 

where  Y is  the  miiltirespnnae  data  matrix 

y = [v,  -y,  y,l 

and  rt  is  the  total  number  of  olxtervations.  Noie  that  in  this  example,  n = 32  and 
p,  = 7.  Using  a FORTR.AN  program  written  by  Conlon  (1992)  (see  Appendix  2.3). 
this  distance  function  is  minimized  and  the  simultaneous  optima  of  all  three  responses 
over  It  are 


The  minimum  of  given  above  is  the  one  obtained  by  using  the  square  root  tran&> 
formation:  to  find  the  minimum  in  terms  of  the  original  scale  of  measurement  of  ps, 
we  simply  si|uare  1.39  and  get  1.83.  Tliese  simultaneous  optima  are  attained  at  the 
location  xi«^  = -1.34  and  x^osi  = —1.21,  tiiat  is,  a holding  period  of  1 week  at  0*C. 
At  this  location  we  have 

l’nrfi7j(x.^)|  = 0.11, 

Vnr(ii,(*„„]]  = 0.30. 


The  design  used  iu  this  example  blocks  orthogonally.  From  Theorem  2.5.2  we  then 


Tlie  second  term  in  tlie  above  equstton  is  equal  to  a zero  matrix  since  ^ 2?=i 
Q.  This  tiieaiis  tliat  'here  is  no  iiirrease  ui  the  predirtioii  \ariaiire  due  to  blorkinq. 
general,  using  Tiieorem  4.5.2,  we  ran  say  tliat 

Vfirlq(j!)|  = Vnr(q,(*l| 

provided  that  m = n^  = . . = To  prove  this,  we  liOte  that  if  ni  s nj  s . . . 
n,  = n‘,  then  the  (|uantity  given  by  (cLi  2;)  " x]  simplifies  to  0. 


CHAPTER 


MULTIRESPONSE  MODELS  WITH  RANDOM  BLOCK  EFFECTS 


III  Clmiiuir  4,  a discussion  of  multircsponse  modeis  with  block  eRecls  was  ^ven 
with  the  assumption  that  the  block  effects  being  considered  are  fixed.  It  is  not  always 
the  case  that  the  block  effects  in  an  experiment  are  Axed.  There  are  many  experimen- 
tal situations  in  which  the  block  effects  should  be  treated  as  random.  For  example,  in 
acbeiiiical  process,  ihe  raw  material  used  may  come  in  batches  in  which  the  quality  of 
the  material  vary  from  batch  to  batch.  In  this  case,  when  we  consider  this  extraneous 
source  of  variation,  the  lilock  effect  (batches)  should  be  treated  as  random.  In  this 
chapter,  we  consider  the  analysis  of  multiresponse  models  with  random  block  effects. 
It  is  important  to  correccty  identify  the  nature  of  the  block  effects  being  considered 
bet:ause  coaductitig  a fixed-effects  analysis  when  the  block  eBeci  is  random  may  lead 
to  incorrect  conclusions.  Tilts  was  sliown  by  Kliuri  (1992)  in  the  single-response  case. 
We  shall  also  consider  interactions  bH.ween  the  fixed  polynomial  effeots  and  the  ran- 
dom block  effects,  that  is.  tlie  block  x treatment  interactions,  in  the  aualysis.  Tlie 
approach  developed  in  Qiapter  4 is  used  to  detennine  the  multiresponse  models  and 
to  obtain  estimators  of  the  model  parametere.  The  results  given  by  Kliuri  (1992, 
1994.  1996b)  are  applied  and  extended  tn  tlie  muiiiresponse  case. 
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lit  i.his5ec[ioitt  wp  follow  tlu*  approncit  riesmbed  In  Chaptpr  4.  Section  4.2  to  estab- 
lish the  niult.itespoiise  nioriel  and  to  deterniine  estimators  for  the  model  parameters. 
Consider  the  same  models  as  in  Seel  ion  3-2,  namely 

V,  = 'Ji>.U  + -V,/3,-Z, «,  + £,,  1 = 1.2 r.  (5.2.1) 

.Assume  that  all  r responses  have  the  same  number  of  observations,  n.  and  that 
6 blocks  are  used  for  earh  response.  Suppose  that  the  block  eHerts  are  random 

such  that  each  6,.  t = 1.2 r,  has  a normal  distribution  with  mean  £(i,)  = 0 

and  E{6fii)  = oa,/»,  i j,  i.j  = 1,2 r,  that  is,  for  the  random  effects  S,  = 

and  6j  = (5^i,  5,j 5,*)',  where  i ^ j,  we  have  ^(5^5^.)  = Ohj, 

II  = 1,2 h and  £(5,,5j„)  = 0,  ii  ?4  u,  u,  c = 1, 2, . . . ,6.  Furthermore,  the  block 

effects  are  assumed  to  be  independent  of  e,  which  also  has  the  normal  distribution 

with  mean  £(e,)  = 0 and  Eie.ef,)  = ir„l„,  i ^ j,  i.j  = 1.2 r.  Note  that  the 

above  models  can  be  rewritten  as 

y.  = .i,,l„^X,0,  + u„  i = l,2,...,r,  (5.2.2) 

where  u,  = Z,5,  -s  a,.  Under  the  assumptions  given  above  we  have 

E(ir.)  = 0,  (5.2.3) 

£(ip,i/')  = o,,/„  -P  at„Z,Z\  = (5.2.4) 

Tlie  muliiresponse  form  of  the  models  given  in  equation  (5.2,2)  Is 

V = iVf S„ -p  X(3 -1- 1-  (3.2.5) 

where  y = y')'.  M = I,  l„,  X = dmg[XuX, X,),  0 = 

i0\.0‘,.  ..0,y.  and  V = (id. id, v',y.  Fs|uation  (5.2.5)  in  turn  can  be  rewritten 


= X"iJ  + j 


(5.2.6) 


whfrp  X“  = (M  . X),  luid  ■d  = {ffu.ffy.  Let  E = (ff|j)  siiri  E<  = (tt»„),  then  from 
equations  (5.2.3)  and  (5.2.4)  we  have 

E[v)  = 0 (5.2.7) 

£(>/*/')  = (E  N /„)  + Z(E,  M /»)Z'  = ♦ (5.2.8) 

wliere  Z = diag{Z],  Zi, Z,) 

To  obtain  estimates  of  the  parameters  in  mode!  (5.2.5)  we  use  the  generaiiaed 
least  squares  (GLS)  approarh  as  desrribed  by  Sriva.stava  and  GileB(1987,  p.l2).  The 
GLS  estimates  are  given  by 

j = (X“'*-'X")-‘X“'*-‘u,  (5.2.9) 

/'do  ^ ^ \ 

Idj  x'*-‘x  j 

Tile  parameter  ^ is  of  interest,  tlierefore.  we  will  concentrate  mainly  on  its  estimate. 
Tlie  GLS  estimate  of  3 from  e<|aaiion  (5.2.10)  can  be  obtained  by  using  the  result 
given  by  Seber(1984,  p.519.  formula  .43.1)  fur  finding  the  inverse  of  a partitioned 
matrix.  .Applyiitg  tide  result,  the  GLS  estimator  is  given  by 

X'»-‘M(Af'$-'M)- 


(5.2.11) 


The  rouditioi]  for  orthoRunal  bltii’kiiift  in  the  imjliiresponse  cAse  And  its  advamaf^es 
were  given  iu  Section  4-n.  This  roiidition  will  be  used  in  the  present  section  to 
determine  its  effect  uii  the  GLS  estimate  of  $.  Tlie  effect  of  orthogonal  blocking  on 
the  parameter  estimates  was  discussed  by  Khuri(1992)  in  the  single-r«tponae  rase.  He 
shovreri  that  if  the  design  blocks  orthogoiiallv,  then  the  GLS  tstimate  of  0 when  the 
block  effects  are  random  Is  et|ual  to  the  GLS  eatiiiiate  of  0 when  the  block  effects  are 
fixed  and  biah  are  etjuai  to  the  OLS  estimate  of  0 (whicli  ignores  the  blork  effect). 
We  now  extend  the  development  to  show  that  the  GLS  estimate  of  0.  when  the  block 
effects  are  random,  is  the  same  ae  the  GLS  estimate  of  0 when  the  block  effects  are 
hxed  in  the  multiresponse  case. 

From  the  results  given  by  Zelliisrl  lBGS),  we  know  that  3 given  in  etjuation  (5-2-11) 
is  the  BLUE  of  0.  If  we  ran  show  that  the  GLS  estimate  of  0 when  the  block  effects 
are  treated  as  fixed  (when  iu  fact  they  are  random)  is  the  same  as  )3.  then  the  former 
estimator  will  have  the  same  optimal  properties  as  0.  This  means  that  when  using 
the  GLS  estimator  in  (4.5.12)  we  are  not  using  an  estimator  that  is  less  efficient  than 
0. 

Suppose  that  the  block  effeeus  are  treated  as  fixed  (when  in  fact  they  are  random). 
Then,  the  variaiice-ruvariance  matrix  given  Iu  eijuation  (S.2.8)  simplifies  to  £ M 
.Assuming  that  the  design  blocks  orthogonally  and  using  Theorem  4.5.1,  the  GLS 
estimate  of  0 can  be  derived  from 

M « IJy. 

0 = jX'V-'X  - X'V-’MlM'V-’Mr‘M'V-‘X]-‘IX'V's,- 


(5.3.1) 


(5.3.2) 


whfre  V = 2 N Equiuioil  (5.3.2)  follows  from  equation  (5-2-11)  by  replacing 
♦ by  V.  To  show  that  3 (GLS  estimator  when  block  effects  are  random)  given  in 
etpiation  (5.2.11)  simplifies  lo  0 (GLS  estimai.or  when  block  effecu  are  fixed)  given 
in  equation  (5.3.2)  under  orthogonal  blockiitg  we  need  to  consider  several  identities 
which  are  established  in  the  following  lemmas'. 

Lf.mma  .7..V,  I Let  A = £4  /&.  then 

= (V  + ZRZ'r'  = V*'  - v'-’ziR-’  + z’v'zr'z'v-'.  (5.3.31 

= \M‘V-'M)-'  + {M'V-'M)-‘M‘V-'ZHR-'-i-Z‘V-'Z- 
Z'V-'‘MtM'V-'M)-'M'V~'Z]-'Z‘V''M{M'V-‘M)-' 

= (M'V-'Mr' -t-C.  (5-3.4) 

where 

e = lM'V-'Mr'M'V-'Z[R-’  -I-  Z’v-Z  - Z'V-'  M{M"V'M)''M'V'Z\'' 
Z'V-'M(M'V-'M)-^  (5.3.5) 

Proof:  Tbe  above  results  follow  directly  from  Seber  (19S4,  p.519,  formula  A3.3). 
Lemma  5.  '/.a  Lj  a d&iijn  Alocits  nrthoqnnaliy.  then 

X‘V'M(M'V-'M)-'M'V-'Z  = X‘V-'Z.  (5-3-6) 

Proof:  Recall  that  M = 1.  -.  1„  and  V = £»/„.  Then  (M'V'M)-'  = 

From  these  results  we  have 


X‘V'MlM‘V-‘Mr'M'V-'Z 


X'V*'M--£M’V-'Z 


_ X‘V-'-(ESc JOV-’Z 


= X'V-'i(E:sJ„)(S-''./,.)Z 

= X'V-'-JZ.  15.3.7) 

wlierp  J = Tr  ts:  J,.  L‘sins«iuation(4.5.9),  hiivp  X'V''^JZ  = X'V'Z.  Notp 
that  «iu»tioii(4.5.9)  is  trup  provided  that  the  rifsign  blodts  orthogonaliy. 


Z‘V-'M{M'V‘'M)-'M'V-'Z  = Z'V-'-JZ.  (5.3.8) 


Proof:  Using  The  results  from  Lemma  5.3.2  we  have, 


Z'V'MtM'V-'Mr'M'V'Z  = Z'V'M-ZM'V'Z 
= Z'V'i(Ex  J,)V-'Z 
= Z’V'ifES;  J„)(E*’ x/,)Z 
= Z'V'-JZ.  (5.3.9) 

Lfinmu  uti)=  -.(fl-'+Z'V'Z)-’  ondEs  (flr'+Z'V'Z-Z'V'iJ’Z)-', 

!>  + £!  = |R-'  ^Z’V-'Z)-'Z'V-'-JZ{R-'  Z‘V-' Z - Z‘V-'ijZ)-'. 

(5.3.10) 

Proof:  The  above  result  foliows  from  an  identity  given  by  Searie(l!W?.  p.l51. 
Exercise  16(e)). 


Ill  Chapter  4,  Theorem  4.5.1.  we  have  derived  an  iniportam  result  on  the  GLS 
estimator  of  0 mider  ortlioKomd  liloi-kiiiB  wiieii  rlie  block  effects  are  considered  to  be 
zero.  In  the  following  ilieoreiii.  we  demonstrate  another  advantage  of  using  designs 
that  block  orthogonally,  if  the  design  used  for  a multiresponse  experiment  witli 
random  block  effects  blocks  orthogonally,  then  'lie  GLS  estimator  of  5 is  not  affected 
by  an  incorrect  analysis  wherein  the  the  block  effects  are  (incorrectly)  considered  as 
fixed  effects.  I.'sing  Lemmas  5.3.1  to  5.3.4  we  .shall  now  prove,  rbat  0 = 0 in  the 
following  theorem. 

Thrnirm  .i..V  I If  the  dcsiyii  klocfct  orlhogontllu.  lAeii  the  GLS  estimate  of  0.  when 
the  block  effects  are  treated  as  fixed,  is  eijual  to  the  GLS  estimate  of  0 when  the  block 
effects  are  random,  that  is.  0 = 0. 

Proof:  Let 

A = .V'*-'X  - (5.3.il) 

8 = (5.3.12) 

then  from  et|uation  (5.2.11 } we  have 

0 = A-'B.  (5.3.13) 

I.'sing  equation  (5.3.3)  in  Lemma  5.3.1  and  equation  (5-3.11),  we  oau  rewrite  A as 
A = A:->-A:1.  (5.3.14) 


a;  = x'*-'x 

ST  X'(V-' - 

= A-i-A,, 


Z'V-'Z)-‘ZV-‘]X 


(5.3.15) 


(5.3.161 


Ax  = X‘V-'X 
A^  = -XV-'Zffi-' *Z'V-Z)-'Z'V-’.Y  (5.3.17) 

a:  = -X'<>-‘M{M'i-M)-'M'i-'X  (.5.3.18) 

= -Y'[V-‘  - V'ZIR-'  + Z'V-'Z)-'Z'V-']Af{M'lV-'  - V'ZlR-'  + 
Z'V-'Z)-'Z‘V-  ]M)-'M'[V-'  -V‘ZlR-‘  + Z'V-'Z)-' 

ZV-'IX.  (,i.3.19) 

From  ^uatioti  (5.3.4),  wp  know  that 


wlirre  C ia  da5n«l  In  equation  (5.3.5).  Tills  implies  that 

A = -X'|V-' - V-'Z(«-' + Z’V-'Z)-'Z'V-‘|A#[(M'V'-'Mr' +C] 


- V-'Z{R-'  -r  Z'V-‘Z)-'Z'V-']X 
= A,  + A.-t-  ..^Aio  (5.3.20) 

Ax  = -X'V-'Af(Af'V-'JVf)-'M'V-'X  (5.3.21) 

A,  = X'V-'MiM'V-'M)-'M'V-'ZiR-^-^Z'V-'Z)-' 

Z'V'X  (5.3.22) 

A^  = -X'V-'MCM'V-'X  (5.3.23) 

A = X'V-'MCM'V-'Z{R-‘ *Z'V-'Z)-'Z'V'X  (5.3.24) 

A,  = X'V-‘Z(R-'-pZ‘V-'Z)-Z'V  M(M'V-'Mr' 

M'V'X 


(5.3.25) 


A.  = -X'V-'ZtR-'-Z‘V'Zr'Z'V‘MIM'V-'Mr'M'V-'Z 

Ift-'  + Z'V'Zr'Z'V-'X  (5.3.26) 

A)  = X'V-‘Z(R-'--Z'V-'Z)-^Z‘V-'MCM'V-'X  (5.3.27) 

Aio  = -X'V-'Z{R-'  *Z‘V-'Z)-'Z'V-'MCM'V-'Z(R-‘  - 

Z'V-'Z)-'Z'V-'X.  (5.3.28) 

To  bp  able  to  rewrite  0 in  riie  desireri  form-  we  apply  Lemmas  5.3.2-5.3.3  to  equa- 
tions (5.3.22 )-( 5.3.28)  to  get  the  following  results: 

A = X'V-'MiM'V-'Mr'M'V''ZIR-'-i-Z’V-'Zr‘ 

Z'V'X 

= X'V-'Z{R-' ~Z"V-'Z)-'Z‘V-''X  (usinp  equa(ion(5.3.6)) 

= -A  (5.3.29) 

A = -X'V-'M{{M'V-'M}-'M‘V'Zl(R''  + Z’V-'Z]-Z'V-'M 


(u.ri'nj  equn(ion(5-3-5)) 

= -X'V-‘Z[{R-'  +Z’V-'Z)  - Z'V-'MiM'V-'M)-'M'V-'Z]-' 


Z‘V-'  MiM‘V-'M)-  M'V-'X  (uatng  equoiton(5.3.6)) 


= -X'V-‘Z\(R-'  -t-  Z'V-'Z]  - Z'V-'M{M'V-'M)-'M'V-'Z]-' 
Z'V"' X (using  ('guaticm(5.3.6)}  (5.3.30) 

= -i.  Z'V-'Z)  - z'v-'-yz]-'z'v-'x 


titstng  egTtati£m(S.3.8))  (5.3.31) 

Using  the  following  re.siilt  used  in  the  aiinplifirstion  of  A.  that  is. 

X'V-'MC  = X'V‘Z\{R-' +Z'V-‘Z)~  Z'V-'MtM'V'Mr  M V-'Z]-' 


Z’V-'MIM'V'M)-'. 


CAQ  simplify  ^ lo  iht*  ronti 


(M'V-'Afr  M'V-'ZiR-  - Z‘V-'Zr'Z’V  X 
= X‘V-'Z[[R--  -r  Z‘V-'Z)  - Z‘V-'-JZ\-'ZV'-JZ(R-'  + Z'V'Z)-' 
Z'V-'X  (t..nnsr9iin(ion(5.3.8)).  (5.3.32) 


From  Kjuaiions  (5.3.25 )-{5.3.28)  wp  ha^'p 


A-,  = X'V-'Z[R-  -i-Z'V-'Zr'Z'V-'X  l,usingegmitwn{i.3.S))  (S.3.33) 
A = -X‘V-'Z(R-'  + Z'V-'Z)-'Z'V-'-JZ{R-''  ■^Z'V'Zr' 

Z'V~'X  (using a^iaf*on(5.3.S))  (5-3-34) 

A = X'V-'ZiR-'  +Z'V-'Zr'Z'V-'M{[M'V-'M)-'M‘V-'Z{R-'  + 
Z‘V-'Z  - Z‘V-‘M(M'V-'M)-'M‘V-'Z]-'Z'V-‘M{M’V-'Mr'} 
M'V~'X  fu.smp fgu<i(.iDn(5.3.5)} 

= X‘V-'ZiR-'  - Z‘V-'Zr'Z‘V-'-JZ\{R-'  - Z'V-'Z)  - 

Z'V-'-JZ\-‘Z'V-'MiM'V-‘Mr'M'V-'X  (using  egtiafitm(5.3.8)) 
= X'V-'ZtR-'  + Z’V~'Z)-'Z'V-'-JZ{{R-'  * Z'V-'Z)  - 

Z'V~‘ —^Z]-' Z'V~‘ X ( using  sgun/t£m(5.3.6(]  (5.3.35) 

Ao  = -X'V-'ZIR-' +Z‘V-‘Z)-'Z’V-'M{{M’V-'M)-'M'V-'Z 
|(H-’  ^Z'V'Z)  - Z‘V-'MlM'V-‘M]-‘M'V-'Z]-'Z'V'M 
-^Z'V-'Z)-' Z'V-'X 

l«sitigsg«n!ion(5.3.5)) 


= -x'v-'ziR-'  + z'v'zr'z'v-'-jz[{R-'  ^ z-v-'Z]  - 


Z'V-' -JZ]-'Z'V-'-JZ(R-  + Z'V'Z)-'Z'V-'X 

{using  f'^nn2ion|5.3.8)].  15.3.3G) 

Cousiderviiuaiioiis  (5.3.31) -(5.3.36)  for  i.lip  terms /4b  to  /tic,.  We  can  further  simplify 
by  looking  at  the  following  groupings  of  tliese  terms: 

A-A  = -X'V-'Z[{R-'  + Z‘V-'Z)-Z’V-'-JZ\-‘ 

[/j  - Z'V-'-JZ{R-'  + Z‘V'Z)-’]Z’V-'X  (5.3.37) 

A-e/4.  = X‘V-'Z{R-' + Z'V-'Z)-'[It- Z'V-'-JZ{R-' *■ 

Z'V-‘Z)-']Z'V-'X  (5.3.38) 

As  + Au,  = X'V-'ZlR-’’  +Z'V-‘Z)-‘Z‘V-'-JZ[{R-'  + Z'V-'Z)~ 


Z'V'X  (5.3.39) 

Tlius  using  equations  (5.3.14)  - (5.3.17),  (5.3.20).  (5.3.21),  (5.3.29),  and  (5.3.37)  - 
(5.3.39),  we  liave 

A — At  + ^ A^  — Ai  + {As  + /4e)  -e  (*4?  A^}  + {As  + pIio) 

= At  -t-  At  +■  [As  -t-  A)  + {At  + A)  + {As  + -4io), 


A = X’V-'X  - X‘V-'MlM'V-'M)-'M'V-'X 


X'V-'Zl(R-'  4.  Z’V-’2fl  - Z'V-'-JZ]-'[I,  - 
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Z'V-'-JZ^R-'  -Z'V-‘Zi-''Z'V-'X 


- X'V-'Z(H"  +Z'V-'Z)-'Z'V‘-J’Z[(R’'  -t-Z'V-'Z)  - 
ZV-'-JZ\-'\h,  -Z'V-'-JZIR-‘  +Z'V-'Z)-‘l 
ZV'X  (3.3.J0) 

Combining  and  grouping  similar  rarms  wp  havp 


- iR-'  + z'v-'zr’z'v-'-j'z((R-'  +z'v-'z)  - z‘v-'~yzj-'i 
;ii,-  z'v-‘-jziR-‘  •i-z'v-'zi-'iz'v-'x.  15.3.41 ) 

Similarly,  wp  can  rpwrirp  B ns  followp 


- X'V-‘ZI~UR-‘  +Z'V-‘Z)  - Z'V-'-J’ZI*'  +(R-‘  -Z’V-'Z)-' 

- iR-'  + Z'V-'Z)-‘Z'V-'ij-Z|(fl-'  +Z'V'Z)  -z'v-'ijz|-'j 

iIt-Z‘V-'^JZ(R-'  ^Z'V-'Zr'jZ’V'i/.  (.5.3.42) 

.Applying  Lemma  5.3.4  to  equation  (5.3.411,  wp  liave 

A = X'V-'X  - X'V-'MiM'V-'M]-'M'V-'X  + 


X'V-‘Zl~E  - -P  + (P  + - 

z'v-'-jziR-  -I-  z'v-'zr'jz'v'x. 


(5.3.431 


SimiLarW, 


B = X'V'y- X'V’-‘MfM’V-'M)-‘M'V-‘y  + 

X'V-'ZI-E  - P + (P  + B)|[7.  - 

Z‘V-'-JZIR-‘  ^Z'V-'Zr']Z'V-'v.  (5.3.44) 


Tlius  from  KiURtions  (5.3.43),  (5.3.44),  and  (5.3.13)  it  follows  that 

0 = \X‘V-^X-X‘V-'M{M'V-'M)-'M'V-'X\-'[X‘V 
X'V'MiM'V'Mr'M'V-'y] 

= 0. 


(5.3.45) 

(5.3.46) 


The  imaraclions  between  tlie  fixed  polynomial  eHecis  and  the  random  effects  Si, 

Sr  are  ignored  in  model  (5.2.1).  Tn  this  section,  these  interaction  effects  are 

added  to  file  model.  Tills  introduces  additional  random  effects  into  the  mndels.  In 
this  case,  the  resulting  process  variance  contains  terms  for  the  block  effects  and  the 
polynomial  interaction  terms.  The  interaction  terms  are  essentially  the  products  of 
corresponding  polynomial  terms  (elements  of  X,)  and  block  effect  terms  (elements 
of  Z,).  Funhermore.  since  the  Z.s  consist  of  ones  and  aeios  only,  the  values  of 
tliese  products  basically  depend  on  tiie  eiemenus  of  X,  (whicli  in  turn  depend  on  the 
settings  used  for  the  input  variablesl.  ,4s  a result,  ilie  process  varianre  will  depend 
on  ilie  settings  of  the  input  variables.  In  addition,  since  ihe  process  variance  contains 
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interaciioii  terms  then  it  follows  that  these  interactions  ran  slfect  the  determination 
of  iiptimuiii  ronditions.  Heiite.  the  addition  of  interaction  terms  between  the  fixed 

polynomial  etfecLs  and  the  random  effects  4|,  4j S,  should  be  considered.  The 

niodei  for  the  ith  response  with  the  added  interaction  terms  is  pven  by 

y,  = .ikl,+X,j3,  + ^,4,  + ET,»A,.  + e„  i = l,2 r,  (5.4.1) 

where  .io,,  X„  0,,  Z„  4„  and  e.  are  as  defined  in  model  (5-2-1),  T,»  is  amatrix  of  order 
tr  X 5 whose  fth  column  is  obtained  by  multiplying  tlie  elements  of  the  fcth  column 

of  X,  with  tile  corresponding  elements  of  the  fth  coiurrm  of  Z,.  k — 1,2 p,, 

f = 1 , 2,  - . . .5,  and  A,t  is  a vector  of  interaction  coefficients  of  order  5x1  associated 
with  the  ail  column  of  X„  k = 1.2 Model  (5.4.1)  can  be  rewritten  as 

y,  = .3o,l„  + .¥,A+’i’,  f = l-2 (5.4.2) 

or  ei|uivalently. 

y,  = f/,2,  + to.  1 = 1,2 r.  (5,4.3) 


U,  = (1„  : X,].  (5.4.4) 

=.  = (db.,/3;)',  (5-4.5) 

to,  = Z,4, +T(A(  + e.,  (5.4.6) 

T;  = [T„  :T„:  ...  :T„J,  (5.4.7) 

A-  = (a;,,a;, a;,,)'.  (5.4.g) 

Since  the  block  effects  are  random,  then  6,  and  A)  are  random  vectors  (i  = 1,2 ,r). 

.Assume  the  following: 


,V(0,  ‘ = 1,2, 


r,  k = 1,2, . 


\ = r*  = l,2 p„  J = 1.2 pj: 

• £(A,kA;,)  = ^ 0 i = j,  i.j  = \.2 r.k^l: 

[ ,=j,  i,J  = l,2 rk  = l. 


hen«,  E(A;a;  | = Z'>  where  T'>  is  a matrix  of  order  p,  x p,  given 
by  E"  = icra).  . j.i.j  = 1.2 r.  k = 1.2..-.?.,  I = 1.2 py.  and 


. a;.  and  e.  are  in 
. a:  and  S,  are  indrp 


lependentfor.j!j,i.j  = l,2...., 
= r;  and 

i^yij  = 1,2 r. 


r„  = Tr.,z,z;  + IT(E'^»/.)T;'+.T.,J„  i54j.i,j  = i. 
r„  = rrr„Z,Z;+T:(E'‘NA)T:'  + o../„  i = l,2 r. 


■y'J.  V 


and  from  (5.4.9). 


E(W)  = Z(£i  - /i)T‘' +(£«/,)  = ? (5.4.12) 

whnre  Z = diap(Zi,  Zpi.  T'  = riiaglTl.Tj, . - . ,T'}.  £«  — and  £a  is 

llip  partitioned  matrix  £i  = [£’').  The  GLS  estimator  of  5 is  given  by 

= = (l/T-’t/r'l/'r-'v  (5.4.13) 

and  its  variance-rovarianre  matrix  is 

l'fir(3)  = (UT-'y)-‘.  (5.4.14) 

The  variance-covariance  matrix  F is  usually  unknown  and  therefore  must  be  es- 
timated in  order  to  obtain  the  solutions  in  (5.4.13)  and  (5.4.14).  A consistent  esti- 
mator of  r can  be  obtained  by  following  the  procedure  described  by  Singh  and  Ullali 
(1974,  pp.192-193).  This  procedure  Is  essentially  a multirespoiise  analogue  of  the 
Hildreth  and  Hourk  (1968)  proredure  for  single-response  models  with  random  coef- 
ficients- The  procedure  involves  the  establishment  of  a regression  relation  between 
the  ordinary  least  squares  residuals  obtained  from  each  fitted  tth  response  model  and 
a linear  combination  of  the  matrices  given  by  U,,  Z.,  T,',  and  In  - Cf*(I/Itf,)“'  U[, 

I = 1.2 r.  which  are  used  to  obtain  these  residuals.  In  this  regression  relation. 

the  paraineiem  of  interest  are  the  varianres  of  the  random  effecla  in  the  model.  We 
illustrate  this  procedure  for  the  multiresponse  rase  by  following  the  steps  outlined  by 
Raj  and  I'llali  (1981.  pp.  241-237);  First,  we  obtain  the  ordinary  least  squares  (OLS) 
residuals  for  eadi  etiuation  in  (5.4.3)  and  write 

ii;,  = Af,i;„  f = l,2 r (5.4.15) 

will. . \f,  = /,  - U,{U;U,)-'U;.  Sintf  M,U,  = O.  then 

M.tti,  = M.(v,  - t/,E,)  = M,p,,  I = 1,2, r. 


(5.4.17) 


I m,(z,z;m,„,  * T„r„M,.«  + . . . * T.,r„  m,,;;,  * . . 
l*T,„r,M,,,.,*M.,.), 


= [mlZ,Z’,m„:m'„T,XjtTn,i:...-.m',,T„T„m,^:...: 


(5.4.21) 


TP  M„U*  j)  is  s matrix  of  order  a’  x [p.p,  + 2)  of  rank 

M,  = l(m;',)' : (m?,)' : . .. : «,)’ : .. . : (rr*)']',  : 

M„  = [(mV,)' ; (mV,)'  : . . . : (m",)'  ; . . . : < = 1 


(5.4.23) 
T- 2 given  by 
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Oncf  all  the  ff„'s  (i,  j = 1.2 r|  arf  obtained.  «•?  can  determine  the  eistimates  of 

the  matrices  Z,  and  £a*  Substituting  these  e.stimated  values  in  et)uation  (5.4.13) 

f = Z(2»  I,)Z'  + T'(Sa  sv  /,)T-'  + (E  kI„).  (5.4.28) 

Note  that  otiier  estimatore  of  T can  be  obtained.  For  example,  GLS  estimation  may 
be  applied  to  equation  (5.4.26).  Other  approaches  suggested  by  Srivastava  and  Giles 
(1987,  pp.  336-337)  include  obiainiug  ir;^  as  a ininimum  norm  quadratic  estimator. 
Using  r in  place  of  F in  equation  (5.4.13),  an  EGLS  estimator  of  = is  given  by 

=.  = {U't''U)-'U‘t''y.  (5.4.39) 

Singh  and  Ullah  (1974.  p.l93)  showed  that  E,  is  consistent  and  that  the  asymp- 
totic distribution  of  nj(S,  - =1  Ls  normal.  Hence  an  approximate  expression  for  the 
variance-covariance  matrix  of  S,  is  given  by 

Var(5,)  = (t/’f''y)-'.  (5.4.30) 

j.5  Mean  Resuoiise  and  Prediction  Variance 

The  mean  of  the  tth  response  at  a point  x = (i|,  Z], . . . ,xt)'  inside  the  experi- 
mental region  92  is  given  by 

q,(x)  = ?;(x|=,  i = l,2 r,  (5.5.1) 

where  c((x)  = (l./I(x)],  /((x)  is  a vector  of  the  same  form  as  a row  of  X,  iu  (5.42), 

but  is  evaluated  at  the  point  x,  i » 1.2 r.  From  (5.5.1)  and  using  the  results 

from  Section  4.4.  the  predicted  value  for  tlie  mean  of  the  ! th  response,  i = 1 , 2 r, 

at  a point  x = (X|,  . . . ,Ck)'  in  9i  is  of  the  form 


(5.5.2) 


where  2^  is  obtained  from  equation  (5.4.29).  The  multiresponse  form  of  equa- 
tion (5.5.2)  is  ^iven  by 

y,(x)  = e'(*l^  (5.5.31 

where  Q{x)  = dini7(5i(®).5.j(il ?,(*)).  The  estimated  predirtion  varianre  ma- 

trix is  approximately  etptal  to 

'■or[i;,(i)|  = g'lxm't''u]-'S{x).  (5.5.4) 

5.6  Diitimizatioii 

To  determine  the  optima  of  the  r responses  of  interest,  we  use  the  multiresponse 
upcimization  procedure  described  in  Section  3.3.  The  simultaneous  optima  of  the  r 
responses  ran  be  determined  by  finding  rnnditions  on  x that  minimize  tlie  distance 
fiinctian  given  liy 

= ily,(*l  -<f>r|ff'(®}H^'f’V)-'ff(x)r’l!/=(®)  -t^D*  (5.G.1) 

where  the  elentents  of  <j>  are  the  individual  optima  of  the  elements  of  y,{x)  over  the 
experimental  region 

The  individual  optima  and  the  ininimum  value  uf  p[yc(x).0|  ton  he  obtained 
by  usiug  the  toniroiled  random  search  pnn'edute  presented  by  Price  (1977).  A For- 
tran program  that  implements  the  procedure  has  been  written  by  Conlon  (1992).  A 
description  of  Price's  procedure  as  implememed  by  Conlon's  Fortran  program  (first 
mentioned  in  Gtaptcr  2)  is  given  in  .Appendix  B,  Section  B.3. 

The  steps  to  be  taken  to  find  the  minimum  of  the  distance  function  given  iu 
equation  (5.6.1 ) are  as  follows: 

I Estimates  of  F 2a.  and  2a  are  obtained  according  to  equation  (5.4.27)  to 
deterniiiie  the  estimate  of  T froin  (5.4.23). 
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of  3 U9  then  obtained  using  (5.4.29). 


3.  The  estimated  responses  are  optimized  individually  to  obtain  the  vector  4 of 
estimated  individual  optima.  The  loiitrolled  random  search  prorediire  given  by 
Price  can  be  used  to  obtain  these  individual  optima. 

4.  The  distance  function  is  minimized  over  S.  .4gain,  Price’s  procedure 

can  Ire  used  lo  ohuin  this  minimum. 

0.7  ri'Weness  to  Threet  \'aiues 

[n  some  pzpenmemal  situations,  it  may  be  desired  to  determine  the  settings  of 
the  independent  variables  that  result  in  a predicted  response  close  to  a target  value 
I.  The  simpirsc  way  to  measure  the  closeness  of  the  predicted  responses  to  the  target 
values  is  to  use  the  proceiiure  described  in  Section  2.3,  that  is,  we  use  a distance 
function.  Let  t,,  i = 1,2,,. ..r,  be  the  target  value  for  the  ith  response  aod  let 
e = We  measure  the  closeness  of  !/,(x|  to  t by  using  the  distance 

function  given  in  ecfuaiion  (5.6.1},  with  replaced  by  t.  Tbe  distance  function  in  this 
case  is  given  by 

= lly.(»)  - tris'(r)(t;'r‘'u)-'g(a)i-'(y,(*)  - t]}i.  (5.7.1) 

Tlius  in  order  lo  achieve  a certain  degree  of  closeness  to  f.  values  of  x should  be 
determined  such  that  the  dlstmice  funrtion  given  in  (5.7.1)  is  minimized. 

j.a  A Numerical  F.vaninle 

For  this  example,  we  use  the  results  of  the  experiment  by  Chen  et  al  (1969) 
mentioned  in  Section  4.G  concerning  the  post-storage  condition  of  freshly  harvested 
Delicious'  apples.  In  this  case,  we  consider  apples  harv’esied  at  the  same  date  from 
four  differwii  urcharris.  The  four  orchards  were  selected  at  random  from  all  the 
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orchards  located  in  the  Mid-Columbia  district  in  the  Pacific  Northwest  of  the  United 
States.  Thus,  orchard'  can  be  considered  as  a random  block  effect.  Recall  that  the 
input  variables  were  A*i  * number  »f  weeks  in  storage  after  harvest  and  -Ye  = storage 
Temperature  (*C)  after  harvest.  Tlie  respoii.se  variables  were  ’ji  = flesh  finmiess  (N), 
j/3  = amount  Ilf  extractable  juice  (nil/lOOg  tliiid  weight),  and  — internal  ethylene 
composition  of  the  apples  (ppm).  Recall  that  a 4 x 4 factorial  design  was  used.  The 
ilesigii  settings  in  the  original  and  coded  variables  and  the  multiresponse  data  are 
given  in  Table  5.1.  Note  that  the  square  root  of  the  values  of  pj  are  given  In  this 
table.  We  use  the  .same  ending  for  X|  and  given  in  Section  4.6.  that  is, 

■Y,  - 2.S 
" 1.12 
X,  - 9.3 
“ 7.70  ■ 

From  Table  5.1  we  observe  that  all  the  responses  have  the  same  blocking  arrangement. 

Here,  as  in  Section  4.0,  we  wish  to  obtain  optimum  conditions  on  the  input  van* 
ahles  storage  time  and  temperature  for  apples  liarvested  from  orchards  located  in  the 
Mid-Coliimhia  district,  that  is,  find  the  simultaneous  optima  of  the  three  responses 
taking  into  account  the  block  effecis  represented  by  the  different  orchards.  To  do 
this,  we  first  set  up  the  multiresponse  model  so  that  EGLS  estimates  can  be  ob- 
lained.  Once  these  estimates  are  detennined,  we  use  the  .simultaneous  optimisation 
procedure  given  in  Chapter  3.  A complete  second-order  model  in  i|  and  is  fitted  to 
each  response  data  under  the  assumption  that  the  full  model  is  adequate  for  all  three 
responses.  .A  st|uare  root  i raiisformation  was  used  to  stabilize  the  residuals  for  y^. 
This  means  that  the  vector  used  in  the  model  and  parameter  estimation  contains 
the  stiuare  roots  of  the  observed  values  for  ip.  The  model  for  each  response  is  given 


y.  = u.,z, 


1.2,3, 


ORCHARD  .Vi 


Orisinal  laput  VnririilKi  Coded  Input  Variablw 


0.45  51.77 
0.45  41.77 
0.45  38.2E 
1.38  53.67 


59.51  73.5  2.74 


-0.64  57 
■0.64  55 
0.45  55 
0.45  52 


05  74  3-7881 


75  71.5  10.57 


40.07  72.5  11.4^ 


19.41  73  13.83 


Table  5.1.  (continued) 


14  74.5 


97  71.S  12 


•■5.1749 

4.5299 

7.0548 


97  69  14.4069 
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where  = [Is,  : X„],  X.  is  a matrix  whose  iiih  row  is  given  by  /'(x„)  = 

« = 1-2 64,  ui,  = Z„S,  + T;a;  + e,,  Z,  = T^  k lu, 

T;  = ;T„i  : T,«  : Trf  ; Tm  : Trt).  and  (k  = 1,2,3, 4.5)  is  a 64  x 4 matrix  whose 
(th  ((  = 1,2, 3,4)  eoiumn  is  obtained  by  muitiplying  the  kth  column  of  X„  with  the 
corresponding  elements  of  the  ftli  column  of  Z„. 

Using  equation  (5,4.27)  we  And  the  estimates  of  the  matrices  S,,  Sa,  and  S. 
Tlieae  are  given  in  ,Appendix  F Using  these  estimates  we  obtain 

Hi  = (50,8250,-4,6323,-7,5035,0,36505,0,9969.-3,3048)' 

Hj  = (72-3761,-0-5481,-1.5790,0,2727,-0,0771,-0.4880)' 

Hs  = (7.5013,0.8648,3.7189,-0.2439,-0.3410,0.5803)’. 

Using  the  above  values,  we  have 

j.,(*)  = 50.8250  - 4.63231, -7.50351}  4- 0.365051?  4- 0.9969i5-3.3048x,x, 
«rt(a)  = 72.3761  - 0.54S1T,  - 1.5790x,4-0.2727i?-0.077l:c? -0.48801, 1, 
j,d(*)  = 7.5013  4- 0.86481,  4- 3.71891} -0.24391? -0.34101?  4- 0.5803iii,. 

It  is  desired  Lu  iiiaxiinize  both  v,  and  p},  and  minimize  ps  over  the  experimental 
region  S.  Here,  fi  is  the  region  given  by  -1.34  < xi  < 1-34-  —1-21  < Xj  < 1-38. 
The  minimum  of  the  distance  function  given  in  equation  (5.6.1)  is  found  by  using  the 
Fortran  program  presented  by  Conlon  (1992)  (see  Appendix  B,  Section  B-3).  The 
simultaneous  optima  obtained  are  as  follows: 

/ 62.7374 
= 74.5787 

V 1.9291 

located  at  i,  = -1.3178  and  ij  = -1.1928,  that  is,  a holding  period  of!  .02  w-«ks  at 


0.12'r. 


CHAPTER  6 


DESIGN  CRITERION  TO  DECREASE  THE  INTEGRATED  MEAN  SQUARED 
ERROR 

li.l  Iniroiiuction 

The  Iasi  three  chr^rlers  have  dealt  mainly  with  multirespoiise  uptimizaiion  and 
aualyses.  Response  surface  methodology  does  not  only  deal  with  problems  of  opti- 
mization and  analyses,  but  also  includes  techniques  for  the  selection  of  designs  for 
which  certain  criteria  are  satisfied-  In  this  chapter,  we  consider  the  problem  of  find- 
ing a design  for  which  a multirespoo-se  version  of  the  integrated  mean  squared  error 
(IMSE)  introduced  by  Box  .uid  Draper  U959]  is  minimized.  To  minimize  the  IMSE 
we  consider  the  reduction  of  the  variances  of  the  predicted  values  of  the  responses 
and  at  the  same  time  seek  to  establish  sums  protection  against  the  possibility  that 
the  fitted  model  is  an  inadequate  representation  of  the  true  surface.  Thus  we  con- 
sider both  variance  and  bias  errors.  Tlie  resulLs  given  by  Box  and  Draper  (1959)  are 
for  single-response  models  only.  Kim  and  Draper  (1994)  extended  this  criterion  to 
the  multirespoiise  case  assuming  that  the  variance  covariance  matrix,  S.  for  the  r 
responses  is  known.  Tlieir  results  were  limited  to  the  rase  of  two  response  variables 
,ind  unly  one  input  variable  of  interest.  In  this  chapter,  we  extend  the  IMSE  criterion 
to  the  multirespoiise  case  witliuut  making  the  assuinplioii  tliat  £ Is  known.  We  apply 
some  of  the  results  given  hy  Kim  and  Draper  (1994)  to  establish  the  inultiresponse 
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IMSE  and  present,  a sei|uemial  promiure  that  works  for  tiny  number  of  responses  and 
input  variaiiles. 

f>  2 \Tndpl  and  \nratiml 

Let  n tie  tlie  total  number  of  experimental  runs  and  r be  the  number  nf  response 
variables  of  interest.  We  assume  tliat  all  the  r responses  depend  on  the  same  set  of 

input  variables,  Jti.ui, Xt.  and  they  I’an  be  represented  by  a polynomial  of  dejtree 

fi  within  a region  of  interest,  R.  For  the  ith  respon.se  the  fitted  model  is 

U,  = X0,  + e„  1 = 1,2,.. ,,r,  (6.2.1) 

where  }(,  is  an  n x 1 vector  of  observations  on  the  ith  response,  X is  ao  n x p matrix 
of  rank  p of  known  functions  of  the  settings  of  the  input  variables,  ,8,  is  apx  1 vector 
of  unknown  ronstniit  parameters,  and  e,  is  tlie  n x 1 vector  uf  errors. 

The  model  for  the  true  ith  response  mean  (i  = 1.2 r)  is  assumed  to  be  .a 

polynomial  of  degree  d with  possibly  some  additional  terms  of  degree  higher  than  d. 
It  ran  be  written  as 

E(V, ) = ,,,  = Xfl.  + Z7„  1 = 1.2... .,r,  (6.2.2) 

w)iere  £(11,)  = rj,  represents  the  expected  lalue  of  y,  under  the  true  model,  Z is  an 
n X ^ matrix  of  known  fuiictiuns  of  the  settings  of  the  input  variables,  and  7,  is  a 
ij  X 1 vector  of  unknown  parameters.  If  the  fitted  model  of  (6.2.1)  is  correct.chen  7, 
in  (6.2.2)  is  equal  to  the  zero  vector. 

The  multiresponse  form  of  the  models  given  in  (G.2.1)  and  (6.2.2)  can  be  written 


£(y]  = 7 = XB  + Zf, 


(6.2.4) 
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rfspeflivpiy,  whfrf  Y = \y^  ■■  y-,  y,],  B = |fl,  : j3j  : . . . : fl,],  T = [7-1  ; 7,  : 

. . 7,j.  and  « = [«i  : «2  - - : e,)-  Thp  rows  of  y are  indppp:itl?ni  willi  a common 

variAiicp-covariance  malnx  £ = (cr,,)  of  ortiar  r x p. 

The  r x 1 vectors  of  preriicteii  values  and  the  rorrespondiiig  true  mean  values  at 
some  point  i inside  the  experiinenial  region,  ft.  are  given  by 


«!(*)  = BrW)'  (6.2.8) 

»?(*)  = (Till*).  rh(®), . . . • hrfa))'.  (6.2.6) 


respectively,  vrhere  n(»)  = /'(*)(3„  t — 1.2 r,  j/ifx)  = /'(ili3„  i = 1-2 r, 

/'(z)  is  a vector  of  the  same  form  as  a row  of  X but  is  evaluated  at  the  point  z, 
and  0,  is  the  ordinary  lea.st  squares  estimator  of  0^  given  by  0,  - (X  X)"'X  y,. 
Note  that  the  ordinary  least-squares  estimator  of  0,  is  the  same  as  its  generaliaed 
least  squares  estimator  (see  Srivastava  and  Giles  (1987.  pp.  17-18)).  Tliis  is  due 
to  our  choice  nf  Htting  all  the  responses  using  polynomials  of  the  same  degree  and 
form.  Fur  the  above  models  and  notation  we  now  turn  our  attention  to  the  integrated 
mean-squared  error  (fMSE). 

G.3  Inieeraieri  Maan-Siiiiared  Error  Criterion 

The  integrated  mean-squared  error  criterion  (IMSE)  is  rnnremed  with  both  the 
variance  error  and  the  bias,  due  to  lack  of  6t.  of  a given  6tted  model.  It  is  described  as 
follows  (see,  for  example.  Khuri  and  Cornell  (1996),  pp.  209-210  for  the  treatment  of 
the  single-response  case):  Suppose  the  Htted  and  true  models  for  the  mean  responses 
ate  given  by  equations  (6.2.3)  and  (6.2.1)  respectively.  Then,  at  each  point  z In  the 
experimental  region  ft,  tiie  .traled  or  weighted  mean-squared  error  (M5E)  of  y{7)  in 
(6.2.5)  Is  dehnrd  as 


.WS£|s)(z)|  = B(|v(z)  -q(z)|'S  ' (v(z)  - q|z)|}. 


(6.3.1) 


If  the  sralcd  A/fiffyll))  is  aversged  over  R.  we  obtain 


H -T|(!r)|'Z‘'(v(«)- t}(a)|}rfz  {6.3.2) 

wtiere  II”'  = jgtix  and  lix  ineaiis  tizi  rttj.. . (fi».  To  be  able  to  compare  different 
designs  which  may  or  may  not  contain  the  same  number  of  points,  tlie  criterion  in 
(6.3.2)  is  put  on  a per  observation'  basis  and  is  expressed  as 

J = rvn^E{[v(®)-rj(i))'E’'[s)(a!l -r|(x)])dx 
= nS!  (r[E{lv(x)  - tj(x)]'E''[jr(x)  - rj(x)|Hcte 
= rdlj^  (r[E-'E((j(x)  - i7(x)||»(l)  - >?(l)r)]dx  (6.3.3) 

where  n is  the  nuntber  of  poinus.  The  difference  y{x)  - t;(x)  can  be  partitioned  into 
i)(x)  - T|(x)  = (j((x)  - E[i;(x)])  + {E[»(x)|  - ri(x)).  (6.3.4) 

This  partitioning  enables  us  to  separate  J into  two  parts,  nauiely, 

/ = nfljf  (r{E''lV|v(x)]}dx 

+ nQ  j^tT{S-'\E[ii(x)]  - rf(x)|[£(v(x)l  - f|(x)l'tdl 
= V + B.  (6.3.5) 

that  is,  the  t)uantity  / can  be  represented  as  the  sunt  of  the  average  variance,  V,  and 
the  average  stpiared  hits,  S,  errors.  Thus  to  minimize  J one  needs  to  minimize  both 
the  variance  and  bias  errors  of  v(x). 

Tile  (juantities  V and  B in  (6.3.5)  can  be  expressed  as  functions  of  the  region 
moments,  design  muments.  and  Z as  follows.  Csing  the  results  given  by  Kim  and 
Draper  (IPtM)  we  can  rewrite  V In  tlie  fnnn 

V'  = nil  J !r(E”'V'nr(y(*)])ifx 
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= n(J  f lr\S-'f'(x){X'Xr‘f[xiS]dx 
= >iSl I ir[ftx){X'xr'ffx)I,\(U 

= nnj^rf{xi[X'xr'f{x)Hx.  (6.3.6) 

Now.  from  Khiiri  anti  Comril  (1996,  pp.  219-220)  wp  have 

V = nriljjrlf'ix^X'xr'fixfiHx 

= nT,Ti[x'xr'\nl^f{xjf{x)dx\] 

= nr(r[(X’X)-'/i.i)  (6’3-7) 

wlierp  Mn  is  » matrix  of  momimua  of  ilm  spherical  region  B pven  by 

=«/^/(*)/'(x)rte.  (6.3.8) 

Note  that  f{malrtx)dx  means  the  corresponding  matrix  of  integrals.  B on  the  other 
hand,  can  be  rewritten  in  the  form 

a = nil  I lr{Z-'t'[Z'xtx'x)-'f{x)  - g{x]][f\x){X'X)-' X' Z - 9'(x)|f  )dx 
= nil  j lr{Z''f[Z'XiX'X)-'fix)fix){X'X]-'X‘Z- 

g{x}f{z){X'Xr'  X'Z  - Z'XiX'x)-'f[x)g'{x)  + ji*)9'(x)]f}dx  (6.3.9) 

where  g'(i)  is  a vertor  of  the  same  form  as  a row  of  Z,  but  is  evaluated  at  the  point 
X.  Using  the  results  given  by  Khuri  and  Cornell  (1996,  p.230),  B 4AD  he  siniplihed 
furtlier  to  tlie  form 

B = nrr(£-'f'(Z^Y(X'X)-'''!  / /(i)/'(x)dx|(X'X)-‘X'.Z 

- |n|^g(x)/'(x)(/xl(X'X)-’x'Z-Z'X(X’x)-’|f2y  /(x)9'(x)dx)  + 
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\{l  l^g[x)^(x)dx])t) 

= nir{T-'t'[z'xix'xr‘^,,{X‘x)-'x'z  - ^',,{x'xr' x' z - 
= nlrlE-'f'Sf'l  (B.3.10) 

where 

’S  = z'x[x'x)-'ii,,{x'x)-'x'z-ii',,{x'x)-'x'z-z'x(x'xr'ti,,*ii^-, 

(G.3.11) 

It,,  and  >(52  ai^  mathees  whose  elements  are  the  moments  of  the  spheriral  region  B 
,t„  = iijj{x]g\xjdx  (6J.12) 

= f2  f^g{x)g'{x}dx.  (64,13) 

If  we  add  and  subtract  lt\,lti,  /t,,  from  "S.  we  get 

^ = (ihi  + ((,*  X)"'X  Z - »sri'/i|2lViil(,y  X)”'X  Z - /inVisl 

(6.3.14) 

B = ntr{£-‘f'if).  (6-3-15) 

Given  the  expressions  for  I'  and  B in  (6.3.7)  and  (6.3.15)  respectively,  we  can  now 

consider  the  minimization  of  the  quantity  J,  given  by 

J = nr  ir|(X'.Y)-Vii]  + n MS-'f'if ). 


(6.3.16) 


shftll  tiow  M»f*k  a klesigti  f<>r  which  J is  tiunimizH,  chat  is,  a J'OplLin&i  desipi. 
y-optiinal  n-puiiii  ilesign  is  one  which  iias  a minimum  value  of  J among  all  ciesigns 
of  Hxeci  size  n.  We  note,  however,  thai  the  clioite  of  the  design  tliai  miuimizes  J 
depends  not  only  on  X but  also  on  the  matrices  f and  E whicli  are  unknown.  To 
solve  this  problem,  we  propose  to  find  an  expression  that  approximates  J.  dues  not 
depend  on  f and  S,  and  whose  iiilnlniization  will  result  in  a decrease  in  the  value  of 
/ To  do  this,  we  need  the  following  result: 

lrlE-‘r'Af]  <e™.(i)tr(E''f'f)  (6.4.1) 

where  emo«(.)  denotes  the  largest  eigenvalue  of  the  matrix  inside  the  parentheses. 
Proof: 

/r(E-'f'iir)  = frlE-if'AfE-i) 

= irl(fT-i)'ilfE-i)] 

= <rf(fE-i)'(i-e„„(i)/,-ee„^(A|/,I(fr-*)) 

= /.r(-(fE->l'[e™„li)/,-il|fE-^)  + 

-■™lA)(fE-*)'(fE-*)} 

< e^(A)fr[(rE-i)'(r£-i)| 

since  c,n,r(A)/,  - i is  pD.sitive  semidefinite.  Hence, 

<r(E-'f'if)  < e™,lAllr(E-‘r'f ). 

To  Hnd  an  expression  that  does  not  depend  on  f and  E.  we  consider  inequality 
(6.4.1 ).  Using  this  Inequality,  we  have 

nr  lr|(X'A:)-'/»„]  + o rr|£-'f'if]  < nr  !r[(X',X)-Viil  + n e^(Allr(E-‘f'f ). 


(6.4.2) 
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Since  (r(2"‘f'f ) us  conetain,  can  juar  consider  the  following  upper  bound  on  J 
go  eu  by 

7 = c-{r(r|(X'X)-Vn]  + c™(^))  ('*•‘'•3) 

where  c*  = n maxil,/r(E*'f'f  I).  Thus,  a.s  en  nltemative  to  minimiaing  J.  we  need 
only  cuiisider  the  ininlmiuitinn  of  the  ((uuitity  given  by 

r^rlrllX'Xr'/inl-fCe-eii)-  (B-'*'*) 

If  the  niimlter  of  design  points,  n.  is  Hxed  beforeiianri,  finding  a design  tliai  mini- 
miies  J‘  ill  (B.4.4),  with  respect  to  the  nk  design  settings  (coordinate  settings)  in  X. 
is  compuucionily  difficult,  specially  for  very  large  values  of  n or  *.  Thus  we  propose 
a sequential  procedure  wherein  the  design  points  can  be  chosen  one  at  a time.  The 
.sieps  taken  to  find  a design  tiiat  inininiites  J'  are  as  follows; 

I Start  with  an  iuitial  design  D,^,  consisting  of  n.  points  for  which  X X is 
nonsingular. 

2.  Obtain  the  liesign  point  at  which  the  minimum  of  J'  given  in  (6.4.4) 

Is  attained.  Tills  minimum  may  be  obtained  by  using  the  procedure  given  by 
Conlon  (1992)  which  applies  the  optimisation  algorithm  of  Price  (1977)  (see 
.Appendbt  B.  Section  B.3). 

3.  Obtain  D^^.,  by  augmenting  D„,  wiili 

4.  Continue  this  process  to  find  Xn.4j<  Xn„43,  ....  until  the  absolute  value  of  the 
difference  between  the  minimum  of  J'  in  tiie  current  iteration  and  the  minimum 
of  J'  in  ilie  previous  previous  iteration  is  less  than  a small  positive  number 

We  now  illustrate  the  use  of  this  procedure  through  an  esample. 


Talilf  6.1.  Initial  Dwign  1 


We  coosidw  an  experiment  with  two  responses,  y\  and  pj,  which  depend  on  four 
input  variables.  x\,xi,xi,  and  i*.  Tlie  experimental  region  R is  given  by  R = (*  : 
I*  < 4.0}.  Let  the  fitted  morieU  be  a reduced  form  of  a second-degree  polynomial 

U,  = ^ A]X]  + ■fii'i.ciia  + iiiiiii3  -e  -t-  dijjija:,  + 

'fiw-eji, -e  + 1 = 1,2. 

Suppose  that  the  true  models  are  third-degree  polynomials  complete  with  pure  cubic 
and  cross-product  terms  in  and  x«.  Thus  tlie  matrix  X consists  of  a col- 

umn of  ones  and  the  columns  corresponding  to  xi,xs,x.i,  x^.  XiZa.  XiXs,  xix^.  xsx.i, 
X;X4,  and  X3Z4.  The  matrix  Z coiittuns  the  columns  correspondiug  to  i’,  x|,  Z3,  x\, 

X|l|,  Iixj,  Ijxf,  Xtx\,  XjxJ,  X.ll?,  XjX^,  X31},  X4lJ,  X,xl.  Xax],  IlXjX.I,  X|XjI,, 

i,Xji4, 11ISI4,  xj,  xj,  X5,  and  xj. 

Tlie  iterative  procedure  is  carried  out  using  two  initial  designs  given  in  Tables  G.l 
and  6.2.  Tlie  points  in  tliese  two  initial  designs  wen-  generated  randomly  using  a 
iiniform(-l,lJ  distribution.  Tlie  values  of  J‘  for  initial  designs  1 and  2 are  20S7.9308 


Tnblf  6.2.  Initial  Dpsigii  2 


0-2  -0.9  0.1 

1.0  -1-0 
0.1  -0.2 


and  126.4029,  respKi.iviHy.  The  augmented  design  points  and  the  corresponding  min- 
imum values  of  J'  are  given  in  Tables  0.3  and  6.4.  From  Table  6.3  we  can  see  that  the 
procedure  has  been  successful  in  reducing  the  value  of  J*.  We  observe  in  Table  6.3  a 
very  large  liet^rea-se  in  these  miuimum  values,  that  is,  from  179.8272  to  0.79S5.  Sim- 
ilariy,  in  Table  6.4  we  can  see  a reduction,  from  12.7517  to  0.7913,  in  the  minimum 
value  of  J' . It  is  clear  front  these  results  that  ths  choice  of  the  initial  design  has  a 
significant  eHect  on  tlie  locatinn  of  the  new  design  points. 

Suppose  that  is  known  and  is  given  in  column  [1]  of  Table  6.5,  so  that 
lr{E"‘rr)  = 63.4994.  We  can  then  compute  the  values  of 7 and  / at  the  augmented 
design  points  for  initial  designs  1 and  2.  These  values  are  given  in  Tables  6.6  and  6.7. 
Using  column  (1),  we  can  see  in  Table  6.6  (baaed  on  the  augmented  design  points  for 
initial  design  1)  a relatively  large  reduction  in  the  value  of  J.  The  value  of  J was 
reduced  from  126-2755  to  3.6097.  fn  Table  6.7,  using  the  augmented  design  points  tor 
initial  design  2,  we  ubaerve  a reduction  of  the  value  of  J from  54.4934  to  3.3782, 

Other  values  of  fE'i  were  used  to  obtain  the  quantities  7 and  J.  The  matrices 
used  are  given  in  columus  (2)  and  (3)  of  Table  6.5.  The  trace  of  using  colunui 
(2)  is  4272.7836  and  using  rolumn(3]  is  equal  to  768362.86.  The  corresponding  values 
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13  0.S80973  •0.982250 


0.968301  0.928627  0.550561 

0.989971  -0.991625  -0.950226 
-0.965061  0.816088  -0.111830 

0.662220  0.807856  -0.883574 

-0.344738  -0.854130  0.829897 

-0.653085  0.658808  -0.896054 

-0.487555  -0-593758  -0.948731 
.376260  0.809924  -0.855614 


22 

23  -0.0 

24  -0.183724 

25  0.859010 

26  -0.177598 

27  -0.584714 

28  -0.599250 

29  -0.921640 

30  -0.720229 


0.765838  179 

0.710693 
0.735662 


0.673417  0-523824 
0.020541  0.940986 
-0.780538  0.843671 


-0.897498 

-0.902929 

0.750269 

0.690796 


32  0.859010  0.673417  0.523824 


12  0.948174  -0.848618  -0.673076 

13  -0501627  -0.418894  -0.024769 

14  -0-998280  -0,976462  0.264560 


15  -0.547447  -0.572505 

16  -0.990384  0.632S73 

17  0.879863  0.545184 

18  -0.687972  -0.676830 

19  -0.758453  0.981396 

20  -0.112402  0.794415 

21  0.058555  -0.879905 

22  0.536969  -0.806770 

23  -0.973773  -0.986753 

24  -0.686483  0.887854 

25  0.610565  0.978229 

26  -0.49237-1  -0  793021 

27  0.132274  U.880884 

28  -0.998471  0.009779 


U.93IUi:^ 

-0.516979 

-0.89792! 

-0.940951 


0577534 


-0.638103 

0.658523 

0.2U291 

0.235140 


sforlnitialDwign^^ 
X.  minof./- 

-0.495048  12.7517 

0.760178  5-7670 

-0.753714  4.3861 

0.688045  2.9126 

0.663063  2.2456 

-0.807814  1.7543 

-0.112935  1.4424 

0.778010  1.2485 

-0.648637  1.1,522 

-0.313328  1.1075 

-0.985303  1.0578 

0.889255  1.0017 

0.988906  0.9297 

-0.122343  0.9045 

-0.901180  0.8366 

0.149660  0.8017 

-0.926510  0.7913 
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0) 

-1  OnOMl 
-0.783273 


-2.005I2G 

0,4708395 

-1.319051 

-0.978062 

0.2925267 

1-1725859 

0.8042371 

•3.309822 

2.0612323 

-0.555628 

0.1327232 

-1.371969 

-0.238016 

0.1920078 

1.512543 

-0.37045 

0.7037692 

2.0788861 


T«hl>  fi  5 fr  i VViluB 


•1 .324941 
1.8725506 
-1..389G12 
-0.576408 
1.7774064 
-1.119594 
-0.847786 
0.2264552 
0.1106658 
-0.079081 
0.59IS148 
-0.517861 
-0.585036 
0.8573888 
-0.081221 


•0.71061 

-0.153044 

-0.270074 

-0.178132 

0.605084 

0.83527 

1.69032 


11.501244 

-11.76126 

16.283895 

-4.7217.58 

4.7994447 

4.3473334 

6.6230253 

5.3894351 

-9.082806 

0.4155652 

20.003241 

-6.667595 

11.626384 

•4.565.369 

-8.157257 


11.420084 

6.401724 

-0.345166 

2.4910006 

6.9547767 

17.654823 

9.0112095 


2.4076498 

•9.050819 

9.5513472 

-8.774898 

2.3695891 

7.0041361 

11.664103 

•4.508447 

-0.73645 

28.239327 

-6.615139 

15.651509 

6.3082402 

-1.983933 

1.8829675 

0.076124 

13.357089 

2.5683417 

-0.66681 

3.976593 

5.2157095 

17.352086 

-0.452393 


0.2311781  -2.322085 

1.8260022  -2.685598 


(3) 

53.659312  58.020001 
59.327008  55.571656 
41.378685  45-721022 
52.420107  48.320406 
60-270855  51.850057 
45.25225  43.730279 
56.315627  57.057544 
54.678035  62.172334 
43.595999  31.777347 
42.48813  42.151998 
46.59218  55.711063 
34.542682  53.394632 
39.762687  35.599949 
40.605219  35.215676 
42-112611  37.970651 
46.430222  41.323058 
56.917429  65.946676 
42.145796  38.992771 
37.789757  48-730907 
56.9n  584  50.340074 
55.957267  61-82756 

56.309243  50.470919 
43.131956  45.467708 
52.978760  33.778899 
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3 

4 

6 


19 

20 
21 


25 

26 

27 

28 


usini;  ID 
7 J 
6840.13  54.4934 
3329.16  40.544.3 
2849.56  19.2708 
1894.62  16.2182 
1581.08  13.1634 
1293.98  8-6432 


460262.44  2793.7929 

224014.45  2101.9448 

191742.56  1336.4019 
127485.88  1028.372 

106388.81  335.3455 

87070.33  342.0965 


1093.70  5.9545  73727.82 

980.91  5.7843  66004.31 


.305.1616 

132.2252 


070.82  5.2860 

991.90  5.1061 

953.87  4.9745 

050.69  4.9843 

919.30  4.4774 

935.92  3.9967 

901.53  3.9355 

929.82  3.7663 

892.04  3.3782 


6535.31  132.5936 
66743.14  115.4728 
04184.44  123.1381 
G3970.85  118.1996 
61858.08  124.3473 
62976.54  125.0198 
60662.56  123.3794 
62566.51  115.3819 
60023.80  116.2613 


12573714.00  ; 
01197.56.11 
5238133.74 
3482732.70 


2378637.02  24828.663 

2014139.34  22906.978 

1803143.72  1 7170537 

1784594.40  17102.626 

1823327.03  16877.340 

1753427.63  16112-291 

1747592.65  16164.554 

1689874.56  16250.986 

1720429.35  15386.913 

16570554.68  15458. 168 

1709228.03  12384.336 

1639764.90  12368.256 


of  7 «id  J are  alao  fyv»n  in  Tallies  6.6  and  6.7.  Using  coluimi  (2),  wp  obssn.-p  from 
Tables  6.6  and  6.7  a large  reduction  in  the  t'alues  of  J,  that  is.  from  0864.8886  to 
115.51994.  and  2793.7929  to  116.2613,  respectively.  The  samp  is  true  for  the  values 
of  J obtained  from  column  i3).  From  Table  6.6.  wr  note  a reduction  from  244360.31 
to  2501.9079,  and  from  Table  6.7.  231847.53  to  12368.256.  Based  on  these  results,  we 
can  conclude  that  the  procedure  has  been  succea-sful  in  reducing  the  value  of  J. 


The  unknown  values  of  E aud  f make  the  problem  of  minimising  J difficult. 
Box  and  Draper  (1959)  overeame  this  dilficiilty,  in  the  single-respnnse  case,  by  first 

case.  Recall  from  (6.3.15)  that  B Ls  ei|ual  to 


S = nfr[E-'f'Af] 


- - \^x‘x)-'x‘z  - ;.r.V„lVnl(x'x)-‘A:'z  - 


« (rlS-’fV»  ' (r[S->fl(A:'X)-'X'Z-  MuV.,IVu 

l(X'X)-'x'Z-,ir,Vvjfl.  (6.6.1) 


X'z  = ^„. 


1 tliB  prior  liistribuHon.  Thp  design  Is  then  chosen  by  minimiiing  the  expected  value 


CHAPTER  7 


CONCLUSIONS  AND  FUTURE  RESEARCH 


Tills  ilissenatioii  lias  been  primarily  rnncenied  with  the  rieveiopmeiit  of  new  tedi- 
nlques  that  can  be  used  to 

1 set  up  a srrips  of  muliiresponse  expenmenu  (design)  that  will  yield  adequate 
and  reliable  predictions  of  the  responses  of  interest  on  the  basis  of  a certain 
hypothesiaed  model. 

2.  lit  the  hypothesised  inultirespoiise  model  using  data  collected  according  to  the 
design  chosen  in  (1 ). 

3.  determine  optimum  operating  conditions  on  the  model's  input  variables  that 
result  in  the  simultaneous  optimization  of  the  predicted  values  of  the  cepouses. 
within  a certain  region  of  interest. 

In  Chapter  6.  a multiresponse  version  of  the  integrated  mean-squared  error  <IMSE) 
criterion  was  used  to  set  up  a design  that  will  yield  adequate  and  reliable  predictions 
of  the  responses  of  interest.  The  quantity  J'  delined  in  (6.4.4)  was  used  to  minimize 
the  fMSE.  Tlie  iteratii'e  procedure  presented  in  Section  6.4  was  demonstrated  to  be 
successful  in  reducing  J'.  hence  J. 

The  importance  of  taking  into  account  block  effects  in  the  fitting  of  multirrsponse 
models  was  shown  in  Chapters  4 and  5.  When  the  block  effects  are  fixed,  the  effect 
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of  hlodiiiiR  oil  tlip  pfpdifiioii  varifuirp  tnairix  is  given  by  the  difftrentr  Viir(^(i|j  - 
l’nr(^,(s))  whkii  is  posititv  ssmidefinite:  where  mid  ^^(*1  'trp  defined  in  Sertion 
4.4.  Under  the  special  case  of  orrliogoiiaJ  hiocking,  it  was  shown  that  the  generalised 
least  stjuares  estimates  are  the  same  as  when  the  block  efferia  (fixed  nr  random)  are 

To  determine  optimum  operating  roiidliions  on  the  input  variables,  a inultire- 
spouse  optimization  procedure  wherein  the  models  used  do  not  necessarily  liave  the 
same  form  and  degree,  was  introduced  in  Chapter  3.  This  is  a generalization  of  the 
lirocedure  given  by  Khuri  and  Conlon  (1981).  The  procedure  can  be  used  with  any 
number  of  responses  and  any  number  of  input  variables. 

7 2 Future  Research 

Multirrsponse  surface  methodology  i.s  still  In  its  early  stages  of  development.  The 
new  techniques  introduced  in  Chapters  3 to  G cover  just  a very  small  portion  of  the 
work  that  needs  to  be  done  in  this  area  Khuri  (1996a)  lias  given  examples  of  the  areas 
that  need  further  development  and  expansinn  in  multiresponse  surface  methodology. 
Hence  in  this  section,  I shall  only  list  down  some  future  researrh  problems  concerning 
the  results  given  in  Chapters  3 to  6: 

I.  Conlon  and  Kburi(1988)  developed  a [software)  program,  called  MR,  that  rao 
be  used  to  implement  their  inultirrsponse  optiinlaation  procedure  (Khuri  and 
Conlon  (1981)).  This  is  not  the  case  for  the  procedure  given  in  Chapter  3.  Tims 
it  would  be  desirable  If  a software  program,  similar  to  MR,  can  be  developed  to 
implement  the  optimization  procedure  given  in  Giapter  3. 

3.  Most  of  tile  results  giveu  in  Chapters  4 and  5 assume  that  the  same  blocks  were 
iisetl  for ..  ; " esponses.  Hence  there  is  a need  to  generalize  these  results  to 
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the  case  where  the  blocks  used  are  not  neressahly  the  same  for  all  the  response 
variables. 

3.  The  iterative  procedure  presented  in  Chapter  6 to  ininimise  the  inie|rated 
mean-squared  error  is  for  the  case  where  all  the  ntodels  for  the  responses  have 
the  same  degree  and  rorni.  Thus  there  is  a need  to  extend  this  procedure  to 
the  case  where  the  models  used  do  not  have  the  .same  degree  and  form.  Further 
investigations  need  to  be  done  on  the  other  possible  approaches  for  the  mini- 
mization of  J,  mentioned  in  Section  6.6.  Tlie  poasibility  of  applying  the  Box  and 
Draper  |1959)  result  which  suggests  that  a:i  all-bias  design  is  almost  as  good  as 
a design  that  minimizes  J (under  certain  conditions)  to  the  multiresponse  case 
must  be  considered.  .A  procedure  based  on  Bayesian  techitiques  must  also  be 
developed. 

4.  rinally,  all  the  results  given  in  Chapters  3 through  6 assume  that  die  responses 
used  are  continuous.  There  is  a need  to  develop  multiresponse  techniques  for 
the  ronstructioii  of  designs  and  determination  uf  optimum  operating  cooditions 
that  apply  to  experimental  situations  in  which  some  of  the  responses  of  interest 
have  discrete  distributions. 


APPENDIX  A 


MAXaiUM  UKELIHOOOD  ESTIMATOR  OF  2 

A disnission  of  the  estimation  of  the  vahanre-ravariance  matrix,  2.  assoriated 
with  thr.  r responses  of  interest,  is  given  in  Chapter  3.  The  estimator  presented  in 
equation  (3.2.10)  is  that  introduced  by  Zellner  (1962).  Another  estimator  that  ran 
be  used  is  the  one  developed  by  Box  and  Draper  (1965).  This  appproach  uses  a 
Bayesian  argument  and  the  assumptiuii  that  the  random  error  terms  in  the  model 
follow  a multivariate  nonnal  distribution. 

To  determine  tills  estimator,  we  use  the  same  set  up  as  in  Chapter  3.  that  is,  there 
are  r responses  and  k independent  variables  of  interest.  The  relationship  between  the 
response  variables  and  the  independent  variables  is  given  by  the  model 

= + « = 1.2 n.  (A.0.1) 

where  z„  = (iruui-s..-  O'  = (flJ.SJ,  is  a vector  of  ;i  unknown 

parameters,  <ui  is  a random  error,  and  is  a response  function  of  known  form. 


Model  (A.0.1 ) ran  be  represented  in  matrix  form  as 

y = f(£>-.e-)-t-«-,  (A.0.2) 

where  Y = jv,  : pj  : . . . : v,)  is  an  n x r matrix  whose  ith  column  is  p,,  the  vector 
i)f  observations  from  the  ith  response,  i = 1, 2, . . . , r,  e‘  = [ei  ; eg  , . . . s»]  is  the 


rorresponding  matrix  of  ranrioiri  errors,  D'  us  a matrix  of  order  nx  k which  represents 
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the  desiRU  matrix  with  rows  equal  to  and  F'tD'.B')  is  a matrix  of 

order  n x r wlinse  ith  column  entisists  of  the  values  of  /‘(in.S')  for  « = 1,2, n, 

I = 1.2, r.  It  is  assumed  that  the  rows  of  e are  independently  and  ideiiticaliy 

distributed  as  N(0.  S). 

Box  and  Draper  showed  that  the  posterior  density  of  &'  is 


when  nomnformative  priors  for  $'  and  S are  used.  Tiie  parameter  estimates  0 are 
rhosen  to  minimixe  |e''«'|-  Tliis  is  oailwi  the  determinant  estimation  criterion. 

Kangand  Bates  (1990,  p,324)  stated  that  tiie  Bayes  estimator  0 whicJi  minimixes 
the  determinant  criterion  is  also  the  maximum  likelihood  estimator  under  the  same 
asssumptions.  Vsinft  this  result,  tiie  maximum  likelihood  estimator  of  £ is  |;iven  by 


M9'|y)xie-'eTi“ 


(A.0.3) 


(A.D.41 


where  i'  is  given  by 


€'  = Y - F'iD',0' 


(A.0.5) 


APPENDIX  B 


SOME  RESULTS  FOR  MULTIRESPONSE  MODELS 

B.l  PRESS  .Uiii  Cu  criiprui  Iiir  Mullirtsuuua!  uioiieis 

Sparks  (1987)  addressed  the  problem  of  deternuning  which  subset  of  tiie  inde> 
pendent  variables  should  be  used  for  preriirting  values  of  the  reponses  fur  seeuiiugl.v 
unrelated  regression  (SLTl)  models.  He  also  derived  the  Cp  and  PRESS  criteria  for 
SL'R  models.  Note  that  in  our  results  %ve  have  used  tlie  estiinaior  5 of  £ defined 
in  equatlou  (3.2.10).  lu  this  case,  we  will  only  consider  the  results  given  by  Sparks 
II9S7)  that  use  this  estimator. 

Consider  the  inultirespoiise  model  given  by 

t,  = C;-r  + 6 (B.I.l) 

where  y,  V.  7,  and  e.  are  as  defined  in  equation  (3.2.4).  Frotn  (3.2.4)  we  recall 
that  r is  the  number  of  responses  and  that  (7  is  a matrbt  of  order  n x q,  where 
q = Sr=i  9.1  ">d  9.  is  defined  in  (3.2.2).  Let  H be  a q x 7 permutation  matrix  such 
that  LfR  = {U,„  : C7«u)-  that  is,  the  columns  of  U are  reordered  sucli  tliat  the  first  p 
I'oiumnsare  the  variables  to  be  included  in  the  subset  SUR  models  p.„  = Cf,„7,„+ein, 
that  is,  p is  the  tiumber  of  variables  in  the  subset  SUR  model  being  considered.  The 
generalized  least  si|imres  estimainr  of  i„  is  denoted  by  7,„(S)  and  it  has  the  same 
form  as  7 In  eqnatiiin  (3.2.5)  except  tiiat  K is  replarrti  by  S and  U is  replaced  by 
Now,  if  £ is  kuown.  the  prediction  based  on  v,„(£)  = tfm 7,^(5!)  is  biased  and 
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SparlcR  (1987,  p.I06,  fHiuation  4.1)  di^n44:  iha  hiaa  as 

S = [/-r  - U„£l7,„ir)l.  (B.1.2) 

inraaiirp  of  tha  ijualit.y  of  prediction  l«!!ed  on  either  tr(!l]  or  y,„(£)  u the  mean 
s(|uare  error  for  prediction  {MSEP),  defined  bw 

,WS£P[v|  = fl(y-y)(v-v)1  (B.1.3) 

= 'V(v)  + Vnr(y)  + E(y  - y)f)y  - ir)'.  (B.1.4) 

Thus  to  compare  the  predictions  baser]  on  y,„(£]  and  y(£).  we  determine  their 
,\fSEP.  Sparks  gave  tlie  following  expressions  for  .VfSEP[y„(E)|  and  ,WSEP[y(£)l. 
that  is, 

-usEPiy„(E)i  = (Ex/„)  + + bb‘ 

:V/S£P(y(E)l  = (Sk  /„)  + 

Based  on  the  above  result*.  Sparks  stated  that  the  prediction  based  on  y,„(E)  will  be 
better  than  the  one  based  ony(S)  if  |.r((E"' to/„)|AfS£P(y,-„lE))-Af5£'P(v(S))]) 
< 0,  that  i.s.  ^(E"'  -v  I„)B  < q - p.  where  q = ^[=1 9n  and  P is  'he  number  of 
columns  of  A question  that  romes  up  to  mmd  is,  how  can  prediction  based  on 
a subset  model  he  better  than  the  full  model'.’  Sparks  (1987,  p.l07)  presented  the 
following  reasona/conditions  for  excluding  variables  from  the  full  model; 

1.  The  predirtioit  bias  is  aero  for  all  the  responses  at  each  observed  data  point, 
that  is. 

\ /,)B  = 0.  (B.1.3) 

2.  The  MSEP  for  y,„(E)  is  smaller  than  the  SfSEP  fory(E)  at  each  and  even,- 
element  of  these  vectors,  that  is. 


fl'(E-‘  --.I„)B<  1. 


(B.1.6) 


3.  The  MSEP  for  is  smaller,  nri  average,  than  the  MSEP  for  y(£), 


BUZ-'  sI„)B<k-p.  (B.1.7) 

Note  timt  all  the  coiulitluns  giveu  above  are  presented  in  terms  of  A = to 

In)B.  But  Z is  unknown  so  we  need  to  Hnd  an  estimate  of  the  above  conditions,  Chat 
bt,  lind  estimates  of  tiie  expressions  given  in  |B.1.3)>(B.1.7),  Sparks  inotimted  the 
estiniation  of  i.hese  expressions  by  hrst  detintiig 

= [s-'  s /„)  - 15-'  /„)[/, ■.  i„). 

.An  asymptotically  best  linear  unbiased  estimator  for  t„„  where  consists  of  the 
eiements  of  7 corresponding  Co  and  consists  of  the  columns  of  U which  are 
not  included  in  the  subset  SUR  model  under  consideration,  is  then  given  by 

Let  e = p - t/7(El.  Sparks  (9R7,  p.l07)  gave  the  following  results:  The  asymptotic 
distribution  of  tlie  statistic  given  by 

ri/T  - II  nr  - nk  — 7im[^iiuti^q^litil'7oi.l 

” c'(S-'b.'/,W(nr-,) 

is  \*(A],  wliere  A = 'v  In)B  is  the  noncenirahty  parameter.  This  means  that 

an  nsymptotirally  unbiased  estimate  of  A is 

F(9-p.nr-<,)-(y-2p).  (B.1.8) 

DeHne 

Cp(5)  = f(?-p,nr-9)-  (9-2p).  (B.1.9) 

Tlien  the  estiiimted  expressions  for  the  coudilion.s  given  in  (B-1.5)-lB-l-7),  using  the 
asyi..,.iu::(aliy  unbiased  estimate  of  A in  (B.1.8).  are 


CplS)  < p. 


(B.1.10) 
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CMS)  < ;i  + l.  (8.1.11) 

CMS)  < k.  (B.1.12) 

rfsperiivfly.  In  our  casp.  wp  ronsider  the  rondilioii  Rivpn  by  Cp(S)  < p bwauFP  wp 
are  more  coiirernecl  with  tlie  prediction  bias.  Thutt  we  want  values  of  CplS)  that  are 
small  and  dose  top. 

Sparks  (1987)  .stated  that  one  critirisui  aRainst  this  Cp-type  criterion  is  that  it 
measurfS  the  fit  of  the  model  at  the  same  data  point*  used  to  estimate  tlie  regres- 
sion roeffirients.  criterion  which  avoids  this  is  the  PRESS  rriterion.  He  gave  the 
following  procedure  for  deiermining  tlie  PRESS  criterion: 

1.  Exclude  the  ith  observation  and  use  the  remaining  observations  to  estimate  the 
regression  coefficients: 

2.  predict  the  ith  response  vector  using  these  estimates; 

3.  compute  the  generalized  Euclidean  distance  of  the  ith  response  vector  from  its 
predicted  value; 

4.  repeal  |3j  over  all  observations  and  sum  these  distances. 

Using  the  procedure  gif^n  above,  let  y(.,)(S)  be  the  predicted  response  vector  at 

the  ith  data  point,  i = 1.2 n,  based  on  n - I observations,  that  is.  with  the  ith 

observation  omitted.  Denote  the  ith  observation  of  tlie  r x I response  vector  by 
Then  PRESS  is  defined  a.s  follows: 

PR£55=  f;[y,„  -M-.itS)!- 

The  advantage  of  this  criterion  is  tliat  it  is  not  based  on  asymptotic  results  unlike 
the  Cp  cKterion.  Sparks  {1981 ) rircuimends  using  tlie  PRESS  criterion  over  that  of 
the  Cp  criterion  when  n ^ q < 25, 


B.2  McEkov'sg 


MrElroy  (1977)  i-niisicIprKi  a measure  of  coodiiess-of-Hi.  in  t.be  form  of  a coef- 
ficient of  decerminaiion.  for  seeuiiiiKly  mirelated  regressions  model-  Srisanava  and 
Giles  (1987.  pages  346-351)  gave  tlie  following  summary  of  the  characteristics  and 
<lerivation  of  McElroy's  R‘ : To  derive  MrElniy'-s  rewrite  the  multiresponse  model 
given  in  (3.2.4),  namely, 

y = U‘r+e. 

ill  the  form 

p = t^  + CfV  + « (B.2.1) 

where  E?  = rfitis(l?i,I?j If,),  l/°  = dio9(l„,  1„, . . . , 1,),  U,  = (17,  : I,).?  is  a 

vector  of  unknown  parameters,  and  7*  is  also  a vector  of  unknown  parameters,  namely, 
the  imercepts.  If  the  multirtsponse  model  is  written  this  way,  then  the  intercept  tenn 
in  each  equation  ran  be  distinguished  from  the  rest  of  the  other  regressors. 

Let  Sbeany  consistent  estimator  of  S.  with  S"'  = A',4,  where  A is  anonsingular 
matrix.  Premulliplying  equation  (B.2.1)  by  (An/„)  and  applying  ordinary  least 
squares  yields  an  estimated  generalized  least  squares  rstimator  of  7,  determined  by 
tbe  choice  of  S,  that  is, 

i.  = ( 7.  ) = W’iS-'  V X /„)!/ 

with  [he  associated  residual  vector 

e = y-U-y,  = y-U^- U‘y°. 

Now,  using  the  ortliogonality  tierweeii  ordinary  least  squares  residue  and  regressors. 

( ) (A'  X /„)(Ax  J„)e  = ( 
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we  gel  U’'(S~'  •;  I,)e  = 0 whirh  implies  that 


U"e  = 0. 


(B.2.2) 


This  means  ihiu  taking  eleviaiioiis  about  the  sample  means  for  the  data  eliminates  the 
intercepts  from  each  equation,  as  in  the  single-response  case.  Let  P = 

.^fter  some  manipulatinn,  we  get  from  model  (B.2.1) 

(/,h.Ply  = (/,s;PlZ57  + e.  (B.2.3) 

Premiiltiplying  eqUBtinn  (B.2,3)  by  (A  N /„)  and  decoiiipusing  the  total  variation 
into  'explained'  variation  and  uiiexpiaiiied'  variation,  we  gel. 

v’lS-’  S P)y  = f U’(S"'  » Pll^  + e'(S'*  ® I„)e  (B.2.4) 

where  the  result  given  by  bs  P)e  = 0 lias  been  used.  A coefficient  of 

determination  for  the  luultiresponse  model  follows  directly  from  equation  (B.2.4); 


^ pm 

y'[S-'<sP)y 


(B.2.5) 


which  ie  equivalent  to 

pj  = 1 _ 

V'(S"'  -c  P)y 


(B-2.6) 


Some  of  the  properties  of  MrElroy’s  H‘  are  as  follows  (see  also  Neudecker(lS91)); 


• Clearly,  0 < < 1 . 


• It  is  invariant  under  changes  of  locatton  and  changes  of  scale  of  y. 

• It  is  related  monutonirally  to  an  appropriate  test  statistic  of  //«  ; 7 = 0 against 
H.:7aS0. 


It  may  he  interpreted  as  the  proportion  of  the  sample  variation  in  (A  >j  J„)y 
explained  l>y  (Aw 
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The  l»*t  result,  given  *bove  shows  the  principal  weakness  of  Me  Elroy  s K‘.  that 
Is.  the  Interpretation  has  only  asymptotic  jiistlHcation  because  5 = (A'A)'’  Is  an 
asymptotic  estimator  of  E.  .Also,  McElroy's  R*  is  derived  on  the  assumption  that 
each  equation  in  the  system  includes  an  Intercept.  If  this  is  not  the  case,  then  the 
ilecompositioii  in  (B.2.4)  which  ensures  liie  K|uivalence  of  the  two  foninilas  in  (B.2.5) 
and  (B.2.6)  no  longer  applies.  Finally,  McElroy's  /?’  nieasure  is  one  for  the  goodness- 
of-fit  of  the  entire  seemingly  unrelated  regressions  model.  If  the  usual  single-response 
coefficients  of  determination  are  (atcuhued  on  the  basis  of  the  generalized  least  squares 
estimates  and  residuals,  the  nsnal  decomposition  of  the  ordinary  least  st|uares  sums 
of  squares  will  not  hold  in  general,  and  ambiguous  measures  will  be  obtained.  It  is 
not  clear  what  meaning  can  be  attached  to  single-rsaponse  goodness-of-fit  measures 
in  tbe  context  of  seeemingly  unrelated  regression  models. 


Conlon  (1992)  presented  a FORTRAN  function.  OPTCRS,  which  globally  mini- 
mises a constrained  function  using  the  Controlled  Random  Search  (CRS)  procedure 
of  Price  (1977).  Tliis  procedure  does  not  require  that  the  function  lia*  derivatives  and 
is  well-suited  for  functions  with  local  minima,  that  is,  a minimum  exists  in  a given 
(local)  region  of  interest.  The  function  OPTCRS  concerns  the  following  minimization 
problem; 


of  finite  volume  in  R"*. 

Conlon  (1992)  dfSaribes  the  procedure  as  follows:  The  controlled  random  search 

algorithm  involves  storing  function  values  (/i /r)  ac  f points  distributed  randoudy 

thruuglioiit  H.  Let  fnt  and  be  die  largest  and  smallest  function  values  in  the 
starve.  Let  and  be  the  rorresponding  points  in  the  parameter  space.  .At 


Minimize  f{9]  subject  to  9 € 8 


I a m-dimensional  rectangle 


each  iteration-  a new  point  $“  is  getierateH  by  refioctinR  a randomly  selected  point 
from  the  storage  through  the  lentroid  of  r;i  other  points  selected  from  the  storage 
without  replacement.  If  /(S')  < /,„„  then  S'  replaces  in  the  storage-  The 
process  is  stopped  when  /,^  and  /,„„  ace  close  together. 

The  steps  of  algorithm  are  giten  below; 

I Load  storage  of  size  / by  generating  points  Sj, . . . ,S(  from  B and  also  storing 

m) m). 

2.  For  each  iteration: 

(a)  generate  a new  point  S'  € B as 


where  is  chosen  randomly  from  Si 6|  and  8 Is  the  centroid  of  m 

points  chosen  without  replaremeni  from  S S|. 

(b)  Check  S'  € 0.  If  not.  repeat  generation  of  S'. 

(c)  If  /(S')  < then  replace  /„„,  with  /(S'),  replace  S„a,  with  S',  and 
determine  new  values  for  and  /,„,„. 

(d)  If  /m„  - /™  < TOL  or  the  number  of  iterations  has  exceeded  a cut-off. 
then  stop. 

3.  Return  /„„  as  the  minimum  value  and  S,„,„  as  the  location  of  the  minimum. 

Conlnn  (1992)  wrote  the  program  using  FORTR.\N  77.  It  consists  of  the  files 
drimr.f,  oplcra.f,  and  opUib-t  and  these  are  .available  from  him  upon  request.  The 
user  needs  only  to  modify  the  FORTR.^N  function  oplfnc  in  driwr./io  suit  his/her 
own  needs.  To  run  the  program  in  the  Department  of  Statistics  romputer  system,  the 
user  must  first  obtain  all  the  needed  files  and  then  compile  the  program  by  typing; 
/77  -0  driver  driver  J uplcra.j  apUii.l.  Then  to  invoke  the  program,  just  type;  driver. 
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For  ilie  first  pxmnpip  given  in  Clmpter  3,  Section  3.4.  the  following  rode  for  driver./ 


double  preciBioo  opters .optfnc ,f  .x.zliai.tflc.tol 

logical  optbad 

exteraal  opcf no , optbad 

dimension  z(3)  ,xliii(2.3}  ,ek(192) 

integer«4  sl,s2 

common  /usesdd/sl,s2 

call  optS9d(&41,437S) 
xlimCl, !)•*! .682d0 

xlimCl,3)— 1.682d0 
zlia(2,l)-1.682d0 
zlis(2,2)-1.682d0 
xUm<2,3)-l.SB2dO 
ni*3 


irmax*100 

itmax-lOOd 

f«optcrs<optfnc,optbDd,z,oz,zlim,np,imaz,itiBax,tol,vk,ier) 
tfrite(6,l0)  f.itmu.ier.  zCl},  x[2),  x(3] 

ID  formate  Function  value  • ’.flO-S.'  in  ',i6.'  iterations. . 
■ ’ ler  ■ ',i3,  ' Located  at  flO.6,  flO.S,  flO.6) 
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logical  fuactioc  optbndtx.az.xlio) 
double  precisioa  z(l] ,xlio(2,l),  radsq 
data  radBq/3.0/ 


if(x(i).lt.zlm(l,i)  -or.  x(i)  .gt.xli«C2,i))  optbad*. false. 

radius  * x(l)»x(l)  ♦ x(2)*x(2)  ♦ x(3)  • x(3) 
if  (radius. gc.radsql  optbud*. false . 


double  precisioa  fuoctioa  optfac(x,nx) 
double  precisioa  x(l),  sigbac 

double  precisioa  betaO.betal ,beta2,beca3,becall.beca22,beca33 
dineasioa  sigbat<2,2) 

data  betaO/-l. 880099294/, betal/O. 097382915/, beta2/0. 009897626/. 
beta3/-0. 009062406/,  betall/0. 102848418/, bat  a22/-0. 001424141/, 
hats33/-0. 003191473/ 

data  sighat(l.l)/l/,8ighat(1.2)/0/.aigbat(2,l)/l/,sigbat(2.2)/0/ 
optfaca  betaO  * betal«x(ll  + beta2»x(2)  + beta3*i(31 

* betall*x(l)*x(l)  ♦ beta22-x(2)'x(2)  ♦ bets33-x(3)-xC3) 


APPENDIX  C 


DESIGN-EXPERT  RESULTS  FOR  CHAPTER  3,  SECTION  3.4 


Numeric  Optimirstion;  FiU  • EXAMPLE!  Run  on  07/28/96  at  16:44:17 


HINIHUH  MAXIMUM  TRANS  MODEL 


1.370  1.920 
20.16  51.59 
16.80  30.60 
29.31  52-70 


None  Red  Quad 
Nona  Red  Quad 
Hone  Rad  Quad 
None  Red  Lin 


FACTOR 


MAX  START  FINISH 

1.000  0.6470  -0.6201 


OPTIMIZATION  PARAMETERS  VLIGHTS 
RESPONSE  SCALE  GOAL  LON  HIGH  1ST  2N0  RESULT  d 
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133 


)-71S6 


-0-6193 


OPTIHIZATIOK  PARAMETERS  WEIGHTS 


..000  -0-7963 
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Cycle  • 10 


Results 


START 

0-421S 

-0-6301 


0-2761 


FIRISB 

-0.6180 


0.9999 


RESPONSE  SCALE 


OPTIHIZATION 


FARAHEHERS 
HIGH  1ST 
2.S  1.00 

SI  1.00 
30  1.00 

60  1.00 


VEICSIS 


1.920  0.S2 
19.13  0.94 
17.57  0.78 
48.98  0.64 


0.7018 


APPENDR  D 


PRO(5FS  OF  THEOREMS  4.5.1  .AND  4.5.2 

Beforp  we  prove  Theorems  4.5.1  luid  4.5.2  we  Hist  recall  some  of  ihe  iiocaiioii  used 
ill  Chapter  4.  We  have  X"  = (M  ; X).  M = /,  S 1„,  2 = diosfZi.Zj, . . .,Z,), 

Z,  » (/„  - ^J„)Z„  J = F = *«!)(/,(*), /si*) /,(*)),  and  G = 

diajlg,(*),3,(«)..--.g,(*))  with  gll*)  = !/!(*) : plj,].  i = 1,2 r. 

D.l  Proof  of  Theorem  4.5.1 

Consider  model  (4.5.10),  that  is.  y = X0  + Zt  + e.  Under  condition  (4.5.5)  of 
orthogonal  blocking  and  by  applying  the  result  given  in  (4.5.9),  namely, 

X’(E-'  •./„)Z  = X'(E-’  ■cJ„)-J’Z.  (D.UI 

the  generalised  lea-st  squares  estimators  0 and  t satisfy  the  normal  equations: 

X'(E-‘  SJ  I„)X0  + X'(S-’  w = X'(Z-'  S /„)y  (D-1.2) 

Z-J"(E-‘sj/„)Xfl+Z’lE-‘--  /„)Zr  = Z'(i;-' .'0/„)y.  (D.1.3) 

!i  is  rlear  that  MM'  = J.  Using  this  result  we  can  now  rewrite  Ih  l.?l  as 

X’(S-'  m7„)M(J-M'ZtI  + X'(I]-'  \/„)Xfl  = X'(S-'  n/„)v.  (D.1.4) 
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If  m now  premiiUiply  ?C|uation  (D.1.3)  by  Ij, 1^.),  wp  get 

M'iT.-'  . /„)^f  - A/'IE-'  /„]X3  = "./jy  (D.1.5) 

since  = M’  and  diagUvlJ, = «'  Using 

the  result 

we  can  rewrite  equation  (D.1.5)  in  the  form 

iVf’(E-'  ■s  h]-JZr  + M'(E-’  y I„)X0  = M'(E-’  « /„)v.  (D.1.6) 

If  we  substitute  MM'  for  J,  we  can  express  equatiou  (D.1.6)  as 

M'(Z-'  \ I„jM(-M‘ZTi  - M'[Z-‘  s /„)A:3  = Af'(E-‘«/„)y.  (D.1.7) 

N'ext.  we  look  at  'hr  cese  where  5 - 0.  If  d s Q,  tlien  the  generalized  lea.st  squares 
estimators  of  and  0 satisly  the  normal  equations: 

Af'(E'‘  N/„)Af3o  + *^’(S"'  = AflE-’s/„)y,  (D.1.8) 

and 

X'(E-'  X/„)M3„  + X'(E-‘  n/„)X3  = X'(E-'  M 7„)y.  (D.1.9) 

By  comparing  etjuation  (D.1.7)  with  equation  (D.1.8)  and  equation  (D.1.4)  with  equa- 
tion (D.1.9),  it  follows  tiiat  0,  which  satisHes  equations  (D.1.4)  and  (D.1.7),  i.s  the 
same  as  0,  wliicb  sattsfias  etjuations  (D.1.8)  and  (D.1.9). 

QJ Proof  Ilf  Theorem  4.5.2 


To  prove  Theorem  4.5.2  we  need  to  coitsider  tlie  following  lemmas: 


fjvnma  I the  motnce.s  X°.  M.  and  F be  defined  as  before,  then 
-M'Zt  *ra  = H'lX-’  lE--  •.  /n)X”]-'X''(E-’  I„)v. 


(D.2.1) 


Proof:  Fruih  i.hf  normal  Hiuations  (D.1.4),  (D.1.7),  (D-1.8)  and  (D.1.91  pvan  in 
the  proof  of  Tlieorem  4.5.1,  we  (set 

M'(E-'  :X)^  I"  ^ = AflE-'  j. 

X'(Z-'  s P,)(M  : X)  (I"  ) = X'(£-'  : X]  ( ) 

The  above  rtsolu  imply  that 

since  the  matrix  [Af  : X|  is  of  full  rulumn  rank.  Thus  if  we  premultlply  the  two  sides 
by  H'  = [/, ; f”),  we  (-et 


/ jo,  + /•,(*! A '1 

1 1 

f 4i'„Zoti4-/;(*)3, 

h 

4i;zoT,+/i(«)/3, 

(D.2.2) 

1 /3br4./;ix)fl,  ) 

1 1 

1 

Note  that  the  left-hand  aide  is  e<|ual  to 


equation  14.5.12),  and  the  right-hand  side  is  K)ual 


ii;,Zn4'i  + \ 

■;l'„Z„Tr  + /',in3,  ) 


It  follows  that 

-M‘Zf  + F‘0  = \ /„)X‘]-'X"'(S-' 


Lrmmfi  Let 

T‘  = '-rt  - -M’Z 

ujActi-  5 = J,  n Ift.  Then 


(D-2.3) 


T'f  = r[Z'(S'’  v/,)Z]-‘Z'|£-'  (u/„ly.  (D.2.4) 

Proof:  Consider  equation  (D.1.3).  If  we  premultiply  both  sides  of  this  equation 
by  nZ'lE-'  yj/niZ]-'.  we  get 

riZ'(£-'  X /„1Z]-' Z N /„)X;3 -I- T't 

= riZ'lE-'  x/„)Z|-'Z'(E-'s/„ly. 

Thus  to  prove  equation  (D.2.4),  we  only  need  to  show  that 

riZ'lE-'  \-/„)Z|-'Z'ij'(E-'  S5/,)X3  = 0.  (D.2.5) 

Tliis  can  be  airotnplished  by  showing  ihat 

riZ'lE-'  \ /,|Z|-'Z'J’'(£-'  N /,)  = O. 


Let  P be  defined 
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wlirrp  n,,  j = 1.2 6.  is  'lis  sise  of  the  jib  Work.  Then,  [Z'(£"'  = 

|E"'  \ = E ^ Next,  note  that  Z'J'  is  a blork-disgonal  inatrbt  of  order 

hr  X nr  given  by  Z' T = /.  s t4 , wliere  Af  is  the  6 x n matrix 


then  E K ?^  = (Aj)  front  these  results  we  obtain  (E(S  = E«"H,  where 

?i  b a 6 X n matrix  consisting  of  elements  ei|ual  to  one-  It  follows  that 

r[z'(E-‘ ..  i„)Z\-'z'j‘{zr' s i„)  = riE  «-k)(e-‘  « /,). 

The  right-hand  side  of  tlie  above  equation  can  abo  be  written  as 

r(E ««)(£"  s/n)  = r(EE-'  w«I,)  =!“(/, ««), 

Using  equation  (D.2.3|  we  tan  rewrite  7“(Tr  WM)  M 

r(Ir  xn)  = jS‘ifr  X 7i)  - X) «). 

To  prove  equation  (D.2,5)  we  only  need  to  show  that  Zt/rXTi). 

Now,  sinre  = l^,  then  i£'|/,  ®«|  = M'.  .•kiso,  ^M'Z(I,»‘H]  = M'  since 

i;,Zo  = (ni,nj,...,n»)  and  i(n,,nj nj|?f  = 1^.  From  these  results,  we  obtain 

riZ’(E-'  '-/„)Z|-'Z'J'(E-’m/„)  = 0 (D.2,8) 

r[z'(s-‘  \ r,)z\-'z'-j'(S-'  Ki„)X0  = o. 
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f.wnmii  r/iecnuarjante  Seltueeri  lAe  fu(m(i(ti!sH'|X"’(E”''v/„)Jr“l"'X  (£  'jo 

I„)y  and  T'[Z'(£-'  • r.,)Z\-' Z'{S-'  . /„(«)  w ifro.  mm  is. 
cov[H'\x•''l^:-'  '--f-jv) 

= o 

(D.2.9) 

Proof;  To  pmvo  e<iuation  (D.2.9),  ii  is  suffirient  to  show  t.liac 

riZ'lE-’  V I„)Z]-'Z‘{^-'  \ /„)X°  = O (D.2.10) 

Lot  N = (ni,Tij, m)'.  Using  flio  result. 

Z’lE"'  >5 1,)M  = £"'  » Z’,1™  = S-'  » N, 

we  get 

r'tz'fE-'  ./,)2|-’z'(£-'  s/,)M  = t’(Es-;f)(e-‘«jv|. 

Now.  since  :FN  = 1,,  we  can  simplify  T'ISNy)(E"'  'tN]  to  the  form  » 

^N)  =r'(/,Nl,)  = T£.  with  T'£  = ^£‘£-^M‘Z£  where  £ is  riehiied  in  Lemma 
D.2.2.  Both  iff  and \M'Z£  are  equal  to  I,  so  that  T£  = O.  that  is, 

r'(Z'(S-‘  N/„)Z]-'Z'(E-'  K/„)M  = O.  (D.2.111 

If  we  pustmultiply  the  two  sides  of  equation  (D.2.8)  by  X and  use  equation(D.I.I), 
we  also  obtain 

T[Z‘{S-'  \ /JZ|-'Z'1E-'  N I„)X  = O.  (D.2.12) 

Etiuation  (D.2.10)  follows  direclly  from  etiuations  (D.2.11)  and  (D.2.12). 

Lrmma  D.I'.i  Recoil  tIuU  T t;no(^.  - mid  f = P,  5>!  1».  TAeti 

/.  £‘(E  X)£  = lE?.i 
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2.  £'(Z  \ r)ZM  = 5E:  <md 

S.  M'Z(En:P)Z'M  = nS. 

Prouf;  Recall  that  E«^  = 1^,)  where  A,,  is  defined  in  equation  (D.2.7).  Note 
that  Z'M  s I,  N.  and  N = (ni.  n? n*)', 

1.  Since  Ex  = (A,,)  and  = o.,Ef=,  tl'^n  £'(Ex^)£  = !£?.,  ^]E. 

2.  £\£\Z)ZM  = £'{SkZ){I,  -N)  = r'(Ecsl»)  = 5E. 

3.  M'Z(ES;F)Z’M  = M'Z(Swli|  = nS. 

Proof  of  Theorem  4.5.2: 

Let  us  write  f)(®}  = CO  as 


The  right-liand  side  of  the  above  ei|uatioii  in  turn  equals  f“3  + \£'r  so  that 


We  apply  the  result  given  in  Lemma  D.2.1  to  equation  (D.2.13)  and  obtain 


fi(*|  = r/3  + ir-r. 


(D.2.13) 


ij(xl  = H'fX-'iE-'  iw/„)X'|-'.Y”'(E-'  + -Af'Z|f.  1D.2.U) 


From  Lemma  D.2.2  we  can 


rewrite  ei|uatioit  (D.2.14)  as 


q|ar)  = H'IX''(E-’*«/,)JC"|-'X''(E-' v/„)y+r[Z'(E-',</„)Z)-'Z'(£-'yo/,)v. 

(D.2.151 
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Using  ei|uatiou  (D-2.15)  and  Lpuima  D.2.3,  ihp  lariaiice^ovariancs  matrix  of  17  is 

= N/niX'j-'H  + riZHE"'  s:/„)Zl-'T, 

r,iri«7(r)l  = lV(^(x)|  + r|Z’(E-'  \ 1„)Z\-'T.  lD.2,161 

Wr  can  rewrits  et|iiiuioti  (D-2-IG)  in  the  fuilowing  form 

rnr[rj(x)j  = V’nrl^„tx)l +T'(E  3>  ^)T 

by  using  equation  (D.2.6).  Consider  the  second  term  on  the  right-hand  side  of  the 
above  equation.  Prom  etjuaiiou  (D.2.3)  given  in  Lemma  D.2.2  we  have 


T(z  \ :f)t  = 

^CiZ  - i£’(£{«.7’)Z’M  ^M'Z(E  N.?)Z'M. 


(D.2.17) 


FinaJly,  we  obtain  tlie  desired  result,  tliat  is. 


rrm(q(x)|  = V«rlq,(x)j-e  ^ 


by  applying  the  results  given  in  Lemma  D.2.4  to  equation  (D,2,17), 


APPENDIX  E 


SAS  RESULTS  FOR  THAPTER  4.  SECTION  4.6 


Tlie  SAS  Syateu 


C(p}  Variablas 


.83953834  0.82847202 


.93063544  0.92320352  4.73352  X1X2 
.84544524  0.S2S8868O  40.02265  X2SQ 
.84023035  0.32311217  42.18285  XISO 


.93654233  0.92714120  4.28666  X2SQ  X1X2 
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). 93132745  0-92115373 

).8461372S  0.82334276 


6.44687  XISQ  X1X2 

11.73599  XISO  X2S0 


0.93723434  0.92516402  6-00000  XISQ  X2S0  X1X2  ••• 


C(p)  Vuiablas  id  Model 


2 0.76214196  0.74673796  27.59794  XI  X2 


80TE:  The  above  variabUs  are  included  in  all  aodeU  to  follov. 


3  0.85492088  0.83937669 
3 0.78304869  0.75980390 
3 0.77093103  0.74638792 


8.69162  X1X2 
24.88691  X2S0 
27.61746  XISO 


4 0.87582760  0.85743169  97T49  X2S0  X1X2 


>.86370994 


>.84351883 


8-71103  XISO  X1X2 


}. 79183775 


XISO  X2SO 


0.75099890  24.90642 


0.88461667  D.B62427S6  6.00000  XISO  X2SQ  X1X2 


N = 32  RegresBion 


for  Depandant  Variebla:  Y3 


CCp) 


Varrablas  in  Modal 


2 0.82314406 


0.81094710 


XI  X2 


NOTE:  Tba  above  variables 


folloe. 


3 0.84646482  0.83001462  6.24679  X1X2 
3 0.83830162  0.82097680  7.86496  XISQ 
3 0.829543SS  0.81128036  9.58032  X2S0 


4 0.86162238  0.84112200  5.26071  XISO  X1X2 

4 0.86286431  0.83106643  6.98607  X2S0  X1X2 


>.84470111  0.82189386 
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5 0.S6802187  0.84264148  6.00000  XISQ  X2SQ  X1X2  ••• 


Re|r«ssioa  Models  (or  Depesdenc  VsriebU:  Y3SQRT 
R-9quare  Adjasted  C(p)  Variables  ia  Modal 


S. 97748  XI  X2 


0.90073792  0.89010269 


0.68833754 

0.88238294 


0.91450272 


0.89614774 


0.92231291 


0.87637371 

0.86978111 


0.90183645 

0.89499966 

0.88076222 


0.90737308 


9.22063  X2S0 
13.37074  XISQ 
16.36360  X1X2 


6.61388  XISO  X2SQ 
8.60675  X2S0  X1X2 
12.75685  XISQ  X1X2 


6.00000  XISQ  X2S0  X1X2 


Note;  Tbe 't 


APPENDIX  F 


ESTIMATES  OF  VARIANCE-COVARIANCE  MATRICES  FOR  SECTION  5.8 


The  followiiiK  i 


' the  matrices  Ej,  E*.  and  E used  in  the  pjtample 


( 0.5654  0.2094  0.2714  A 
0.2094  0.0000  0.0597  , 

0.27U  0.0597  0.0000  j 


( 8.2743  0.9060  -0.9928  \ 

E=  0.9060  1.8794  - 0.0037  , 

-0.9928  -0.0037  1.5333  ^ 


. dtaff(0,0.2574.0,0.0). 


0827  -0 
1321  -0 
0525  -0 
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( 0.0590  -0.0465  0.1250  -0.1280  -0.0468  \ 

0.0543  0.07.56  0.1125  -0.0S49  0.0263 

0.0176  0.0210  0.0441  -0.0928  0.0192  , 

0.0109  0.0004  -0.0542  - 0.0091  -0.0011 

0.0263  - 0.0102  0.0770  - 0.0727  0.0672  / 


( 0.0827  0.1321  0.0525 

-0.0663  -0.0193  -0.0634 
0.1165  0.0672  -0.0691 

-0.1022  -0.02.33  -0.0143 
-0.1238  - 0.0654  0.0138 


0.0136  0.1282 

0.0877  -0.0742 
0.0241  0.0714 

-0-1062  -0.0687 
-0.0353  -0.1049 


= rKnr;(0.1373,O.OI84.0.( 


.0). 


-0.0773 


0.0308 

0.0384 


-0.0262 

0.0015 


-0.0152 


0.0237  \ 
0.0269 
■0.0328  , 

0.0384 
0.0406  J 


rfm9(0,0,0,0,0). 
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